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I .  INTRODUCTION 

One  concept  which  could  potentially  increase  the  firepower  of  the 
"fire-and-forget"  class  of  helicopter  borne  missile  systems  would  be  to 
acquire  and  hand  off  multiple  targets  from  a  precision  pointing  and  track¬ 
ing  system  (PTS)  to  several  missile  seekers  simultaneously,  or  almost  so, 
in  a  short  period  of  time.  A  typical  over-all  fire  control  configuration 
is  shown  in  Figure  1-1.  The  pointing  and  tracking  system  typically  con¬ 
sists  of  an  optics  train,  line  of  sight  (LOS)  stabilization  system,  for¬ 
ward  looking  infra-red  (FUR)  imaging  system,  manual  and  autotrack  system, 
laser  range  finder  and  associated  electronics.  An  imaging  missile  seeker 
could  be  an  infrared  type.  It  is  assumed  that  during  preflight  checkout 
or  during  the  actual  flight,  the  lines  of  sight  of  all  the  missile  seekers 
are  aligned  with  the  line  of  sight  of  the  PTS.  However,  due  to  gyro  drift, 
boresighting  inaccuracies,  vehicle  vibration  and  flexure,  etc.,  the  seek¬ 
ers  will  not  remain  boresighted  with  the  PTS.  Since  the  PTS  has  a  larger 
field  of  view  (FOV )  in  both  axes  as  compared  to  missile  seekers,  it  is 
expected  that  the  FOV  of  all  the  missile  seekers  will  be  located  within 
the  FOV  of  the  PTS.  The  multiple  target  problem  then  becomes  that  of 
locating  targets  and  missile  seeker  aim  points  within  the  PTS  field  of 
view,  deciding  which  target  Is  to  be  assigned  to  each  missile,  generating 
error  signals  to  the  torquers  In  order  to  slew  the  missile  LOS  such  that 
Its  assigned  target  is  In  the  center  of  Its  FOV,  and  Initiating  automatic 
seeker  tracking. 
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Helicopter  missile  fire  control  system. 
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In  general,  the  task  of  locating  a  given  smaller  image  within  a  larger 
image  is  known  as  "image  registration".  The  smaller  image  is  referred  to 
as  the  window  or  the  reference  and  the  larger  image  is  called  the  search 
area.  With  the  above  notation,  in  the  multiple  target  problem,  the  image 
obtained  from  the  PTS  sensor  is  the  search  area  and  the  image  obtained 
from  each  missile  seeker  is  a  window.  Therefore,  there  is  one  search 
area  and  more  than  one,  say  n,  windows.  It  is  assumed  that  all  the  n  win¬ 
dows  are  completely  located  within  the  search  area.  Now  the  problem  of 
multiple  image  registration  can  be  defined  as  that  of  finding  n  subimages 
of  the  search  area  which  best  match  the  n  windows.  Even  though  very  lit¬ 
tle  attention  has  been  given  to  date  to  the  problem  of  multiple  image 
registration,  a  considerable  amount  of  work  has  been  done  in  the  area  of 
single  image  registration.  Several  interesting  problems  such  as  map 
matching,  cloud  motion  tracking,  ship  and  aircraft  identification  are 
solved  through  digital  image  registration. 

In  Chapter  II,  the  problem  of  single  image  registration  is  precisely 
defined  and  varuous  existing  methods  of  accomplishing  digital  image  re¬ 
gistration  are  described.  The  inherent  problem  associated  with  regis¬ 
tration  algorithms  is  their  high  computational  cost.  An  algorithm  which 
is  computationally  efficient  for  single  image  registration  may  not  be 
efficient  for  multiple  image  registration.  A  detailed  comparison  of  the 
important  multiple  image  registration  methods  based  on  the  number  of 
arithmetic  operations  for  software  implementation  and  the  complexity  of 
hardware  for  real-time  implementation  is  presented  in  Chapter  III.  New 
methods  of  accomplishing  multiple  image  registration  which  are  computa¬ 
tionally  more  efficient  than  the  most  conmonly  used  template  matching 


techniques  (correlation  and  sequential  similarity  detection  algorithm) 
are  described  in  Chapter  IV.  Conclusions  and  recommendations  are  given 
in  Chapter  V. 


This  page  intentionally  left  blank. 
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II.  EXISTING  METHODS  FOR  DIGITAL  IMAGE  REGISTRATION 

The  problem  of  locating  a  given  Image  within  a  larger  Image  uses 
techniques  and  algorithms  fundamental  to  the  disciplines  of  Image  pro¬ 
cessing  and  pattern  recognition.  Template  matching  methods  such  as  "cor¬ 
relation"  and  "sequential  similarity  detection  algorithms"  are  widely 
used  for  the  determination  of  local  similarity  between  two  images  [1]  - 
[7].  The  inherent  problem  associated  with  the  above  two  methods  or  any 
image  registration  method  is  high  computational  time.  Several  schemes 
such  as  "two-stage  template  matching"  and  "course-fine  template  matching" 
have  been  proposed  to  speed-up  template  matching  methods  [8],  [S3 .  The 
"method  of  invariant  moments"  which  is  widely  used  in  classifying  an  un¬ 
known  pattern  as  one  of  several  known  patterns  can  also  be  used  to  ac¬ 
complish  digital  image  registration  [10]  -  [14].  Various  methods  of 
accomplishing  single  digital  image  registration  are  described  in  this 
chapter. 


Problem  of  Digital  Image  Registration 
Let  two  images,  S  the  search  area  and  W  the  window,  be  defined  as 
shown  in  Figure  2-1.  S  is  a  MxN  array  of  digital  picture  elements  (pix¬ 
els)  which  may  assume  one  of  G  possible  levels  on  the  gray  scale,  i.e., 

0<S(1,j)<G-l  (2-1) 

for  1  <  i  <  M  and  1  <  j  <  N 
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W  Is  a  KxL  (K  <  M  and  L  <  N)  array  of  pixels  having  the  same  gray 
scale  range. 

0  <  W(*,m)  <  G-l  (2-2) 

for  1  <  t  <  K  and  1  <.  m  <  L 

Let  j  denote  each  unique  KxL  subimage  of  S  whose  upper  left 
corner  coordinates  are  (i,j).  Then  (1,j)  Is  also  called  the  reference 
point  of  subimage  S.  ^  and  the  (M-K+1)(H-L+1)  reference  points  corre- 

»  »  J 

spondlng  to  the  (M-K+1 )(N-L+1 )  possible  subimages  of  S  are  called  allow¬ 
able  reference  points. 

S^U.m)  *  S(i+*-1  ,  j+m-1)  (2-3) 

for  1  <_  z  <_  K  and  1  <_  m  L 

1  <  1  <  M-K+1  and  1  <  j  ^  N-L+l 

When  S  and  W  do  not  differ  in  pixel  resolution  and  rotation  (or 
have  been  preprocessed  to  equalize  the  pixel  spatial  resolution),  digital 
Image  registration  is  a  search  over  the  allowed  range  of  reference  points 
to  find  the  subimage  .*  which  best  matches  the  window  W.  Existing 
techniques  for  registering  an  Image  within  a  larger  image  and  schemes  to 
speed-up  these  techniques  are  presented  in  this  chapter. 

Correlation 

When  two  images  do  not  differ  in  pixel  resolution  and  rotation, 
the  method  most  widely  used  for  image  registration  is  cross-correlation 
[1]  -  [6].  The  elements  of  the  unnormalized  cross-correlation  surface, 
R(1»j)»  defined  to  be 
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R(1.J)  «  2  2  W(i,m)  S..  Ai,m)  (2-4) 

fi\  m*l  1 

for  1  <  1  <  M-K+l  ,  1  <  j  <  *:-L+l 

In  most  cases,  coordinates  of  the  maximum  value  of  correlation 
surface  Indicate  Image  registration.  However,  since  R(1,j)  Is  a  cross 
correlation.  It  Is  possible  that  the  maximum  value  of  correlation  sur¬ 
face  does  not  Indicate  true  Image  registration.  This  Is  Illustrated 
below.  Consider  an  Ideal  case  where  W  exactly  matches  subimage  Sj*  j*. 
Then 

K  L  , 

R(1*J*)  *  I  2  WZU.m)  (2-5) 

1*1  m*l 

A  A 

Consider  a  nonmatching  reference  point  (1,j)  where 


S;  •  max  W(t,m)  *  W 

(t,m) 

for 

A  A 

Correlation  value,  R(1,j),  Is  given  by 


(2-6) 


R(1.J) 


^  I  i, 


(2-7) 


It  Is  easy  to  show  that 

R(U)  >  R(1*.j*)  (2-8) 

Therefore,  even  In  the  Ideal  case  a  search  for  a  maximum  over  the  cor¬ 
relation  surface  does  not  necessarily  yield  true  registration.  In  order 
that  the  maximum  value  of  the  correlation  function  Indicate  true  image 
registration  W  and  j  must  be  normalized.  Elements  of  the  normalized 
cross-correlation  surface  are  defined  to  be 


17 


R(1.j) 


ll  W(i,m)  S1j(£,m) 


for 


c£  £  v^d.iB)]1*  c £  y  s2,  ,(»,»)]’ 

1  <  1  <  M-K+1  ,  1  <  j  <  N-L+l 


(2-9) 


This  obviously  Involves  more  computation  than  the  unnormallzed 
method  given  by  (2-4).  In  spite  of  Its  high  computational  cost,  corre¬ 
lation  Is  widely  used  In  image  registration  for  the  following  reasons: 

1.  Correlation  appears  to  be  a  natural  solution  for  the 
mean- square-error  criterion. 

2.  Digital  hardware  and  analog  optical  devices  implement 
correlation  easily. 

In  general,  the  amount  of  computation  associated  with  any  similar¬ 
ity  detection  method  Is  proportional  to  the  number  of  pixels  In  the  win¬ 
dow  and  the  number  of  pixels  In  the  allowable  search  area.  In  the 
correlation  method  each  of  the  KL  pixels  In  W  is  compared  with  the  cor¬ 
responding  pixel  In  Sj  j  to  compute  R(1,j).  Since  the  correlation  func¬ 
tion  has  (M-K+l ) (N-L+l )  elements  a  total  of  KL(M-K+1 ) (N-L+l )  pairs  of 
pixels  are  compared.  Thus  the  total  computation  time  Is  roughly  propor¬ 
tional  to  KL (M-K+l ) (N-L+l ).  The  approximate  number  of  arithmetic  opera¬ 
tions  required  to  compute  the  normalized  correlation  surface  is  derived 
In  Chapter  III,  Several  modified  versions  of  the  standard  correlation 
algorithm  exist  and  each  version  has  Its  own  advantages  and  disadvantages. 
Vector  correlation,  feature  matching  correlation  and  hybrid  correlation 
algorithms  are  briefly  described  on  the  next  pages 
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Vector  Correlation  Algorithm 

The  standard  correlation  algorithm,  as  a  measure  of  similarity, 
computes  the  cross  correlation  surface  between  the  window  Vi  and  the 
search  area  S  by  using  only  intensity  levels  for  pixels  of  W  and  S.  The 
vector  correlation  method  computes  the  cross  correlation  surface  between 
W  and  S  based  on  gradient  as  well  as  gray  scale  values  of  pixels  [15]. 
Let  S(1,j)  and  G5(i,j}  be  the  pixel  and  gradient  values  of  the  (1,j)th 
pixel  of  the  search  area.  Similarly,  W(i,m)  and  GW{i,m)  are  the  pixel 
and  gradient  values  of  the  (t,m)th  pixel  of  the  window.  Now,  a  two- 
dimensional  vector  consisting  of  intensity  and  gradient  values  can  be 
associated  with  each  pixel  of  S  and  W.  Let  VS(i,j)  be  the  vector  asso¬ 
ciated  with  the  (1,j)th  pixel  of  S,  i.e. , 


VS(i,j) 


Sd.J) 

SS(1,j). 


(2-10) 


for  1  <  i  <  M  and  1  <_  j  <_  N 

w 

Similarly,  V  (i,m)  denotes  the  vector  associated  with  the  U,m)th  pixel 
of  W. 


for 

Let 


VMU,m) 


W(£,m) 

■  GW(i,m). 


1  <  i  <  K  and  1  <  m  <  L 


V*j(i,m)  •  VS(1+t-l,  j+m-1 ) 
1  <  i  <  K,  l<m<L 


(2-11) 


(2-12) 


for 
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Elements  of  the  unnormallzed  vector  correlation  surface  of  S  and  W  are 


defined  to  be 


R(1J)  -  t  I  (VWU,m))T  (vf  Ai, m))  (2-13) 

i-1  m*l  1,J 

for  M-K+l ,  1  <  j  <  N-L+l 

Elements  of  the  normalized  vector  correlation  surface  are  defined  to  be 

£  £  (vwu,»i))T  <v? ,(«.»)) 

KM)  •  -fj? - - -  •  (2-14) 

CZ  I  (Vw(i,m))T  (Vw(t,in))]’s 
el  m*l 


1 _ 


for  1  <  1  <  M-K+l ,  1  <  j  _<  N-L+l 
It  can  be  easily  shown  that 


R(1,J) 


y[WU,m)S.  .(i.m)  +  GW(*tm)  GS.  Ai, m) 
_ j _ _ 

CS  £  tfu..)  *  G«2U,m)]]“ 

m*l 

_ 1 _ 

C2  ICS?  .(t.m)  +GS?  iU,m)]]‘s 
z-1  m-1 


(2-15) 


for  1  <  1  <  M-K+l ,  1  l  j  1  N-L+l 

Since  the  vector  correlation  method  Is  based  on  more  Independent  Informa¬ 
tion  than  the  standard  correlation  method.  It  Is  expected  to  have  a  better 
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performance.  However,  the  vector  correlation  method  Involves  consider¬ 


ably  more  computation  than  the  standard  correlation  method. 

Feature  Matching  and  Hybrid 
Correlation  Algorithms 

Since  Images  W  and  S  are  obtained  from  two  different  sensors,  a 
certain  amount  of  preprocessing  is  necessary  before  computation  of  the 
correlation  surface.  In  addition  to  the  pixel  spatial  resolution  equal¬ 
ization  preprocessing,  if  the  two  sensors  differ  in  dc  gain  and  bias, 
each  image  can  be  preprocessed  such  that  its  mean  pixel  value  is  zero 
and  standard  deviation  of  pixel  values  is  unity.  This  process  is  called 
Intensity  level  normalization.  Normalized  images  can  then  be  correlated. 
There  are  two  variations  of  the  standard  correlation  algorithm  suggested 
by  the  Rand  Corporation  [16].  Based  on  the  preprocessing  technique  used, 
the  algorithm  is  called  either  the  feature  matching  correlation  algorithm 
or  the  hybrid  correlation  algorithm. 

Feature  Matching  Correlation  Algorithm.  In  this  method,  the 
window  (reference)  Is  segmented  Into  homogeneous  regions.  A  homogeneous 
region  Is  defined  as  a  set  of  spatially  connected  pixels  whose  pixel 
values  remain  almost  constant  over  the  region.  Each  of  the  homogeneous 
regions  Is  then  preprocessed  separately  based  on  its  characteristic. 

For  example,  intensity  level  normalization  of  a  homogeneous  region  is 
accomplished  by  subtracting  the  mean  pixel  value  of  the  region  from  each 
pixel  and  by  normalizing  with  respect  to  the  variance  of  pixel  values. 
Similarly,  the  search  area  Is  also  segmented  Into  homogeneous  regions 
and  each  homogeneous  region  Is  preprocessed  separately.  The  preprocessed 
window  and  search  areas  are  then  correlated  using  a  standard  correlation 
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algorithm.  Simulation  results  reported  In  Reference  [16]  Indicate  that 
this  method  yields  a  sharper  correlation  peak  as  compared  to  the  stan¬ 
dard  correlation  method.  This  may  be  due  to  the  enhancement  of  high 
frequency  content  of  each  homogeneous  region.  The  feature  matching 
correlation  also  compensates  for  contrast  reversals  between  the  cor¬ 
responding  homogeneous  regions  of  the  window  and  the  search  area.  How¬ 
ever,  depending  on  the  scene  and  sensor  resolution  it  may  not  be  possible 
to  decompose  all  scenes  into  homogeneous  regions.  Such  scenes  are  call¬ 
ed  non- homogeneous  scenes.  Even  if  a  scene  is  composed  of  homogeneous 
regions,  it  is  extremely  difficult  to  accomplish  segmentation  in  real 
time. 

Hybrid  Correlation  Algorithm.  In  hybrid  correlation  only  the 
window  is  segmented  into  homogeneous  regions.  Each  subimage  of  the 
search  area  is  assumed  to  be  the  matching  subimage  and  is  segmented 
identically  as  the  window.  The  correlation  between  the  window  and  the 
subimage  is  computed  by  matching  each  homogeneous  region  of  the  window 
with  its  corresponding  region  in  the  subimage  and  by  combining  the  par¬ 
tial  results  additively.  Simulation  results  in  Reference  [16]  show  that 
this  method  is  better  than  the  standard  correlation  algorithm,  but  not 
as  good  as  the  feature  matching  correlation.  However,  the  hybrid  cor¬ 
relation  algorithm  has  the  advantage  of  not  segmenting  the  search  area 
and  thus  requires  less  computation  as  compared  to  the  feature  matching 
correlation. 
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The  Fast  Fourier  Transform  Method  of  Computing 
Correlation  Function1 

The  convolution  theorem  of  fourier  analysis  states  that  convo¬ 
lution  in  the  time  or  space  domain  is  equivalent  to  multiplication  in 
the  temporal  or  spatial  frequency  domain.  Since  correlation  is  a  form 
of  convolution,  an  alternate  method  of  computing  the  correlation  function 
thus  exists  [17],  [18].  Let  X  and  Y  be  two  images  of  the  same  size. 

Then  the  cross  correlation  between  X  and  Y  is  given  by 

R(i,j)  «  IFFT{X(U,V)  Y*(U,V)>  (2-16) 

where, 

X.(U,V)  is  the  discrete  fourier  transform  of  X 
Y*(U,V)  is  the  complex  conjugate  of  the  discrete  fourier 
transform  of  Y 

IFFT  signifies  the  inverse  fourier  transform  operation 
The  size  of  the  correlation  surface  is  the  same  as  that  of  X  or  Y. 
However,  R(0,0)  is  the  only  valid  element  and  other  elements  are  ignored. 

Since  W  is  a  KxL  array  and  S  is  a  MxN  array.  Equation  (2-16) 
cannot  be  directly  used  for  the  computation  of  the  correlation  function. 
This  problem  can  be  solved  by  padding  W  with  zeros  as  described  below. 

0  1  S(i,j)  <  6-1  ,  for  0  <  i  <  M-l  ,  0  <  j  <  N-l  (2-17) 

and  0  <  W(z,m)  <_  G-l  ,  for  0  <_  *  <_  K-l  ,  0  <  m  <_  L-l  (2-18) 

Construct  a  new  image  W1  of  size  MxN  by  padding  W  with  zeros. 

!W(z,m)  ,  for  0  <  Z  <  K-l  ,0<m<  L-l 

(2-19) 

0  ,  for  M  >  i  >  K  or  N  >  m  >  L  . 
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Let  S(u,V)  and  W^u.V)  be  the  two  dimensional  DFT's  of  S  and  W. 

Now, 

R(1,j)  -  IFFTCS(U.V)  •  W|(U,V)]  (2-20) 

for  0  <  i  £  M-1  ,  0  <  j  <  N-l 

R(i.j)  is  valid  for  i  =  0,1 ,2,...,M-K  and  j  =  0,1 ,2,. . . ,N-L  and  other 
values  are  ignored. 

To  compute  the  DFT  or  IDFT  of  an  array  of  size  MxN,  MN  log2MN 
complex  multiplications  and  complex  additions  must  be  performed.  MN 
complex  multiplications  are  needed  to  multiply  the  DFT’s  of  S  and  W-j. 
Since  the  FFT  method  yields  the  unnormalized  correlation  surface,  it  must 
be  normalized.  For  large  M  and  N,  the  FFT  method  of  computing  the  cor¬ 
relation  surface  requires  fewer  calculations  as  compared  to  the  direct 
approach.  However,  this  method  requires  an  additional  memory  of  4  MN 
real  words  which  may  not  be  feasible  for  large  values  of  M  and  N.  It 
is  difficult  to  implement  this  method  in  real  time  due  to  hardware  limi¬ 
tations. 


Sequential  Similarity  Detection  Algorithm 
In  this  algorithm,  a  search  over  each  of  the  (M-K+l )(N-L+1 )  refer¬ 
ence  points  is  performed  as  in  correlation.  However,  the  criterion  for 
similarity  at  reference  points  is  significantly  different  from  that  of 
the  correlation  method.  The  unnormalized  error  e'(i,j,z,m)  and  the 
normalized  error  e(i,j,z,m)  between  the  pixel  W(t,m)  and  its  corres- 

sponding  pixel  in  S.  .  are  defined  as 
1  »  J 

e'(i,j,M)  =  |Si  jd.m)  -  W(i,m)|  (2-21) 

e(i,j,z,m)  =  |S,  ,U,m)  -  S,  ,  -  W(i,m)  +  W| 

1  * J  1 


(2-22) 


where 


and 


i  K  L 

4-  Y.  X  W(i,m) 

“■  £=1  m=l 


j  =  kl  X  X  ^ 
lfJ  KL  m»l  1,J 


(2-23) 


(2-24) 


The  correlation  method  yields  the  correlation  surface  as  a  measure 
of  similarity,  while  the  sequential  similarity  method  computes  the  error 
surface  as  a  measure  of  dissimilarity.  The  normalized  error,  E(i,j),  as¬ 
sociated  with  reference  point  (i.j)  is  defined  as 


K  L  , 

E(i,j)  *  t  £  e(i,j,£,m)  (2-25) 

i*l  m=l 

for  1  <  i  <  M-K+l ,  1  <  j  <  N-L+l 

Now  the  problem  of  digital  image  registration  reduces  to  the  problem  of 

finding  S.*  .*  such  that 
■  t  J 

E(i*,j*)  <  E(i,j)  for  1  <  i  <  M-K+1  and  i  t  i*  (2-26) 

1  <  j  <  N-L+l  and  j  ^  j* 


In  general,  computation  of  error  is  simpler  than  computation  of 
correlation  since  addition  takes  less  time  than  multiplication  and  is 
easier  to  implement.  There  are  methods  such  as  Constant  Threshold  Se¬ 
quential  Similarity  Detection  Algorithm"  and  "Monotonic  Increasing 
Threshold  Sequence  Algorithms"  suggested  by  Barnea  and  Silverman  which 
further  reduce  the  number  of  additions  [7].  These  methods  are  based  on 
some  kind  of  guess  work  or  statistical  assumptions  and  cannot  be  gener¬ 
alized.  The  number  of  arithmetic  operations  required  to  implement  the 
SSDA  method  in  its  entirety  is  computed  in  Chapter  III. 
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Moments  Method 


Given  a  two  dimensional  continuous  function  f(x,y),  the  moment 
nipq  of  order  (p+q)  is  defined  by  the  relationship 


xpyq  f(x,y)  dx  dy 


(2-27) 


for  p,q  =  0,1, 2, 3,... 


A  uniqueness  theorem  states  that  if  f(x,y)  is  piecewise  continuous  and 
has  non-zero  values  only  in  a  finite  region  of  the  x-y  plane,  then  the 
moments  of  all  order  exist  and  the  moment  sequence,  {nip^},  is  uniquely 
determined  by  f(x,y)  and  conversely,  {n^}  uniquely  determines  f(x,y) 
[19]. 

The  central  moments  can  be  expressed  as 


for 

where 


pq 


00  00 

*  J  J  (x-x)p  (y-y)q  f (x,y)  dx  dy 


>oo  —oo 


P»q  =  o,i ,2, . . . 

m  mm 

10  -  -21 
m00  00 


(2-28) 


For  a  digital  image  these  moments  are  given  by 

wDQ  =  Z  I  (x-x)p  (y-y)q  f(x,y)  (2-29) 

x  y 

The  normalized  central  moments  are  defined  as 


pq 


i 


for  p,q  =  2,3,4... 


u00 


(2-30) 
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Hu  has  derived  a  set  of  seven  invariant  moments  from  second  and 


third  order  central  moments  [10].  They  are  given  by 

41  *  n20  +  n02  (2~31 ) 
*2  *  ^n20  ”  n02^  ^  ^11  (2-32) 

*3  *  (n3Q  -  3r*]2^  (3n2^  ”  n03^  (2-33) 

*4  *  ^n30  +  n12^  +  ^n21  +  n03^  (2-34) 

*5  *  ^n30  "  3n12^n30  +  n12^n30  +  n12^  "  (2-35) 

2  2 

3 (n21  +  n03^  ^  +  ^3n21  ‘  n03^n21  +  n03^3^n30  +  n12^  " 

(n21  +  n03^ 

*6  =  (n20  '  n02)C(n30  +  r,l2)2  ‘  (n21  +  n03)2]  +  (2"36) 

4nll (n30  +  n12)(n21  +  n03J 

*7  =  (3^12  "  n30^n30  +  n12^n30  +  n12^  '  (2-37) 

2  2 

3(n2]  +  nQ3)  ]  +  (3n]2  -  n03)(n21  +  nQ3)[3(n30  +  n12)  - 

^n21  +  n03^2^ 


This  set  of  moments  has  been  shown  invariant  to  translation,  rotation, 
reflection  and  scale  change  [10]. 

The  method  of  moments  is  used  for  automatic  classification  of  an 
unknown  pattern  as  one  of  several  known  patterns.  Let  a1 ,  a2,  ...,  ak 
be  k  known  patterns  whose  invariant  moments  are  also  known.  Let  x  be 
a  given  unknown  pattern  to  be  classified  as  one  of  the  k  known  patterns. 
Now,  the  moments  method  consists  of  computing  the  invariant  moments  for 
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x  and  comparing  them  with  the  Invariant  moments  of  a-j,  ....  a^. 
x  Is  classified  as  a^  if  the  invariant  moments  of  x  best  matches,  ac¬ 
cording  to  some  prespecified  rule,  the  invariant  moments  of  a^.  The 
moments  method  is  widely  used  in  applications  such  as  visual  and  digi¬ 
tal  pattern  recognition,  recognition  of  two  dimensional  patterns  with 
linear  distortion,  aircraft  and  ship  identification  and  character  recog¬ 
nition  [10]  -  [14].  It  is  not  widely  used  for  image  registration  be¬ 
cause  of  the  enormous  computation  involved.  However,  this  method  looks 
promising  for  the  multiple  image  registration  problem  for  the  following 
reasons. 

1.  If  the  images  do  not  differ  in  rotation  and  resolution, 
it  is  not  necessary  to  compute  the  normalized  central 
moments  or  the  invariant  moments.  The  moment  sequence 
Oripq}  determines  f(x,y)  uniquely  and  can  be  used  to 
characterize  each  subimage  of  S. 

2.  It  may  be  possible  to  obtain  good  registration  results 
using  very  few  moments.  This  trade-off  should  be  in¬ 
vestigated  in  a  follow-on  program  using  simulation  of 
real  digitized  scenes. 

3.  If  the  correlation  or  SSDA  method  is  used  for  multiple 
image  registration,  the  correlation  surface  or  error 
surface  has  to  be  determined  for  each  of  the  windows 
separately.  In  other  words,  computation  is  directly 
proportional  to  the  number  of  windows.  For  the  moments 
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method,  however,  once  the  moments  are  computed  for  all 
subimages  and  windows,  matching  Is  accomplished  with 
negligible  computation. 

As  a  result,  the  moments  method  may  prove  economical  for  the 
multiple  image  registration  problem. 

Hough  Transformation  for  Digital  Image 
Registration 

One  possible  way  of  accomplishing  digital  image  registration  is 
by  detecting  the  predominant  lines  in  the  window  and  the  search  areas 
and  using  that  information  to  find  the  subimage  of  S  which  best  matches 
the  window.  Consider  an  image  consisting  of  a  number  of  discrete  white 
points  (edge  points)  on  a  black  background.  The  problem  is  to  detect 
the  groups  of  col  inear  or  almost  colinear  edge  points  (white  points). 

Of  course,  the  problem  can  be  solved  by  testing  the  lines  formed  by  all 
pairs  of  edge  points.  However,  computation  required  for  n  points  is  ap- 
proximately  proportional  to  n  and  may  be  prohibitive  for  large  n. 

In  1962,  Hough  proposed  an  ingenious  method  of  detecting  lines  in 
binary  images  [20],  [21].  He  replaced  the  original  problem  of  finding 
colinear  points  by  a  mathematically  equivalent  problem  of  finding  concur¬ 
rent  lines.  A  straight  line  is  given  by  the  equation 

y  *  mx  +  c  (2-38) 

where  m  is  the  slope  and  c  is  the  intercept.  Any  line  can  be  uniquely 
identified  by  its  slope  and  intercept.  Therefore,  a  line  in  the  x-y 
plane  maps  into  a  point  in  the  m-c  parameter  plane  and  vice  versa.  Sim¬ 
ilarly,  a  point  in  the  x-y  plane  maps  into  a  line  in  the  parameter  plane. 
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This  can  be  seen  by  letting  x  and  y  be  constants  and  solving  for  m  as 
a  function  of  c  in  Equation  (2-38).  The  result  is 

m  =  —j^c  +  y/x  3  ac  +  b  (2-39) 

where  a  and  b  are  constants.  In  the  m  vs.  c  plane,  (2-39)  is  the  equa¬ 
tion  of  a  straight  line.  Thus,  n  points  on  a  straight  line  in  x-y  plane 
are  transformed  to  n  lines  which  intersect  at  a  common  point  in  the  m-c 
plane.  Therefore,  the  problem  of  finding  col  inear  points  in  the  x-y 
plane  is  equivalent  to  that  of  finding  concurrent  lines  in  m-c  plane. 

Slope  and  intercept  both  being  unbounded  complicate  the  applica¬ 
tion  of  the  Hough  transformation  for  line  detection.  In  order  to  over¬ 
come  this  problem,  Duda  and  Hart  suggested  the  use  of  an  angle-radius 
rather  than  si ope- intercept  parameter  plane  [22].  A  straight  line  can 
be  uniquely  specified  by  the  angle  9  of  its  normal  and  its  algebraic 
distance  p  from  the  origin  as  shown  in  Figure  2-2.  The  line  can  now  be 
represented  as 

x  cos  9  +  y  sin  9  =  p  (2-40) 

where  9  is  bounded  and  takes  on  values  between  0  and  2n  and  p  is  less 
than  or  equal  to  R,  where  R  depends  on  the  size  of  the  image.  From 
Equation  (2-40)  the  following  properties  of  the  Hough  transformation  can 
be  easily  verified. 

1.  A  point  in  the  image  plane  (x-y  plane)  corresponds  to 
a  sinusoidal  curve  in  the  parameter  plane.  This  can 
be  seen  by  letting  x  and  y  be  constants  in  Equation 
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(2-40)  and  solving  for  p  In  terms  of  e.  The  result 
is  p  ■  K  sin  (e  +  0)  where  K  *  +  y^  and  0  = 

arctan  x/y. 

2.  A  point  In  the  parameter  plane  corresponds  to  a  line 
In  the  image  plane.  This  can  be  seen  by  letting  p 
and  e  be  constants  In  Equation  (2-40).  The  result 
Is  y  *  ax  +  b  where  b  *  p/slne  and  a  =  -cose/sine. 

3.  Points  lying  on  a  straight  line  In  the  image  plane 
map  Into  sinusoidal  curves  in  the  parameter  plane 
each  passing  through  a  common  point. 

4.  Points  lying  on  the  same  sinusoidal  curve  in  the 
parameter  plane  correspond  to  a  family  of  lines  through 
one  point  in  the  image  plane. 

Therefore,  if  all  the  edge  points  are  mapped  to  the  parameter 
plane,  the  problem  of  finding  col  inear  edge  points  in  the  x-y  image 
plane  becomes  that  of  finding  concurrent  sinusoidal  curves  in  the  para¬ 
meter  plane.  The  point  of  intersection  of  these  sinusoidal  curves 
uniquely  identifies  the  straight  line  edge  in  the  image  plane. 

Detection  of  Lines  in  Digital  Images 

Let  F  be  a  digital  image  of  size  KxL  whose  pixels  can  assume  one 
of  G  possible  levels  on  the  gray  scale.  Let  GF  be  the  gradient  image  of 
F  which  can  be  computed  using  any  of  the  known  edge  detection  algorithms 
(e.g.,  Sobel  edge  detector  or  Roberts  cross  operator).  GF  is  then  trans 
formed  into  a  binary  image  by  setting  all  pixels  with  gradient  values 
greater  than  a  predetermined  threshold  to  one  (edge  points)  and  all  re¬ 
maining  pixels  to  zero  (non-edge  points).  Let  n  be  the  number  of  edge 
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points  In  the  binary  Image.  Suppose  all  edge  points  are  mapped  Into 
their  corresponding  sinusoidal  curves  In  the  parameter  plane.  In  gen¬ 
eral,  these  n  curves  will  intersect  at  points  corresponding  to 

possible  lines.  Exactly  col  inear  subsets  of  edge  points  can  be 
found,  in  principle,  by  finding  coincident  points  in  the  parameter  plane. 
Unfortunately,  this  method  is  exhaustive  and  computation  grows  quadrat- 
Ically  with  the  number  of  edge  points. 

When  It  Is  not  necessary  to  determine  the  lines  exactly,  following 
Hough's  basic  proposal,  the  p-8  plane  can  be  quantized  into  a  qu^druled 
grid  on  the  basis  of  an  acceptable  error  in  p  and  a.  The  quantization 
Is  confined  to  the  region  0  <_  e  £  2n  and  0  <_  p  <_  R,  where  R  depends  on 
the  size  of  the  image.  Assume  that  an  accumulator  is  placed  in  each  cell 
of  the  grid.  For  each  edge  point  (x^,  y.j),  the  sinusoidal  curve  given 
by  Equation  (2-40)  is  entered  in  the  grid  by  incrementing  the  count  in 
each  accumulator  along  the  curve.  When  all  edge  points  are  mapped,  each 
accumulator  contains  the  number  of  curves  through  It.  A  count  of  k  in 
accumulator  cell  (a^.pj)  means  that  precisely  k  edge  points  lie  (to  with¬ 
in  the  grid  quantization  error)  along  the  line  whose  normal  co-ordinates 
are  8^  and  p^.  However,  the  exact  location  of  these  k  edge  points  on 
this  line  in  the  x-y  plane  is  not  known  (i.e.,  it  is  not  known  if  the 
edge  points  are  adjacent  or  widely  separated).  To  determine  this,  some 
sort  of  connectivity  test  must  be  used. 

Digital  Image  Registration 

Step  1:  Transform  the  window  and  the  search  area  to  binary  images  as 
described  previously. 
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Step  2:  Divide  the  e-p  parameter  plane  Into  a  quadruled  grid  on  the 
basis  of  acceptable  error  In  e  and  p. 

Step  3:  Map  each  edge  point  in  the  reference  into  Its  sinusoidal  curve 
In  e-p  plane.  Even  if  the  acceptable  error  in  e  and  p  is  mod¬ 
erate  (5  degrees  and  3  pixels),  the  parameter  plane  will  have 
more  than  KL  cells.  As  a  result  any  attempt  to  use  all  the 
Information  in  the  e-p  plane  will  increase  computation.  In 
order  to  accomplish  data  reduction  only  the  predominant  lines 
(cells  with  relatively  high  count)  are  retained  as  features  of 
the  window  and  the  remaining  information  is  ignored. 

Step  4:  Determine  the  predominant  lines  present  in  each  subimage  of  the 
search  area  by  repeating  the  procedure  outlined  in  Step  3. 

Step  5:  Match  the  line  features  of  the  window  with  those  of  each  sub- 
image  according  to  some  predetermined  criterion.  One  way  of 
doing  this  in  the  e-p  plane  is  given  in  Reference  [23]. 

It  Is  felt  that  Hough's  method  of  image  matching  is  not  very 
sensitive  to  slight  geometric  distortion  and  rotation  and  is  also  insen¬ 
sitive  to  small  differences  in  pixel  resolution  of  the  window  and  the 
search  area.  However,  this  method  requires  a  highly  reliable  algorithm 
to  quantize  gradient  images  to  two  levels.  In  the  parameter  plane,  cell 
(o^.Pj)  represents  the  line  in  the  picture  plane  whose  normal  co-ordinates 
are  and  p^  and  the  corresponding  count  gives  the  total  number  of  edge 
points  on  that  line.  The  exact  location  of  edge  points  on  the  line  is 
not  known.  Therefore,  in  order  to  determine  a  true  line  ( i . e. ,  points 
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are  adjacent  on  a  line)  and  prevent  Isolated  edge  points  from  affecting 
the  result.  Hough's  method  may  have  to  be  coupled  with  some  kind  of 
connectivity  test. 

Schemes  to  Speedup  Template  Matching 

"Template  matching"  Is  the  comnon  terminology  used  for  correlation 
and  sequential  similarity  detection  algorithms.  In  template  matching, 
each  of  the  KL  pixels  In  W  Is  compared  with  Its  corresponding  pixel  In 
Sj  j  to  compute  the  measure  of  similarity  (correlation)  or  dissimilarity 
(sequential  similarity  detection  algorithm).  Therefore,  the  total  amount 
of  computation  Is  roughly  proportional  to  the  product  of  the  number  of 
pixels  In  the  reference  and  the  number  of  allowable  reference  points  In 
the  search  area.  Since  there  are  KL  pixels  in  W  and  (M-K+l )(N-L+1 )  al¬ 
lowable  reference  points  in  S,  computation  Is  proportional  to  KL (M-K+l ) 
(N-L+l).  Three  popular  schemes  of  accomplishing  savings  in  computation 
and  speeding  up  template  matching  are  presented  in  this  section.  All 
methods  accomplish  savings  in  computation  by  reducing  the  total  number 
of  pixel  pair  comparisons.  These  methods  are  computationally  more  ef¬ 
ficient  In  terms  of  the  number  of  arithmetic  operations  required  (soft¬ 
ware  implementation)  but  may  not  enjoy  any  advantage  in  real  time  imple¬ 
mentation  using  special  purpose  hardware. 

Two-Stage  Template  Matching 

The  two-stage  template  matching  technique,  suggested  by  Rosenfeld 
and  Vanderbrug,  searches  for  j*  which  best  matches  W  in  two  stages 
[8]. 
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Stage  1:  In  the  first  stage,  some  subimage  W  of  size  pxq  from 
W  and  its  corresponding  subimage  Si  s  of  the  same  size  in  S<  .  are 
matched  to  compute  a  measure  of  similarity  (correlation)  or  dissimilar¬ 
ity  (SSDA)  between  W  and  S.  .  for  i  *  1 ,2,. . . ,M-K+1  and  j  *  1,2,... 

•  *  J 

N-L+l .  W  and  S'.  .  are  shown  in  Figure  2-3.  The  net  effect  of  step  one 

»  t  J 

is  to  find  the  (M-K+l)  by  (N-L+l)  correlation  surface  with  a  reduced 
reference  array  size. 

Stage  2:  In  the  second  stage,  all  reference  points,  (i,j)  for 

which  the  measure  of  similarity  is  less  than  (correlation)  or  the  measure 

of  dissimilarity  is  greater  than  (SSDA)  a  predetermined  threshold  T  are 

discarded  as  non-match  points.  At  the  remaining  reference  points,  the 

window  W  is  matched  with  S.  .  in  its  entirety.  The  method  of  finding 

■  *J 

Sj*  j*  is  the  same  as  before.  The  computational  savings  in  this  two- 
stage  recognition  procedure  results  from  not  having  to  match  the  entire 
template  at  each  reference  point.  Savings  in  computation  depend  on  the 
size  of  W  and  threshold  T. 

Rosenfeld  has  suggested  a  method  to  determine  optimal  values  for 
the  size  of  W  and  T  for  a  given  W  and  S.  His  analytical  model  is  based 
on  many  simplifying  assumptions  which  are  rarely  true  for  typical  mil¬ 
itary  type  images.  Improper  selection  of  W  and  T  can  lead  to  the  pos¬ 
sibility  of  discarding  the  true  registration  point  in  the  first  stage 
itself. 

Course-fine  Template  Matching 

Course-fine  template  matching  is  also  a  two-stage  matching  algo¬ 
rithm  [9].  For  the  first  stage,  the  spatial  resolution  of  both  image 


Figure  2-3.  Layout  for  two  stage  template  matching. 
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arrays  is  reduced  by  replacing  each  pxq  subarray  by  its  average.  This 

yields  a  reference  array,  W,  of  size  K/p  x  L/q  and  a  search  area.  S', 

of  size  M/p  x  N/q.  The  amount  of  computation  required  to  find  the  course 

correlation  surface  for  W*  and  S'  is  proportional  to  —  (---+  1 ) 

'  pq  P  P 

N  L 

(-  -  —  +  1).  At  each  of  the  correlation  surface  peaks  for  which  the 
correlation  value  is  greater  than  T,  the  original  reference  VI  is  cor¬ 
related  with  the  original  search  array.  The  largest  value  of  this  cor¬ 
relation  computation  is  treated  as  the  registration  point.  Savings  in 
computation  depend  on  p,  q  and  T.  This  technique  can  also  be  used  with 
the  SSOA  method. 


Hierarchical  Search  Method 

This  technique  is  a  generalization  of  the  "course-fine  template 
matching"  scheme.  In  this  method  the  search  for  S-*  ,*  is  done  in  n- 
resolution  levels  [24],  [25].  From  a  given  window  W  and  the  search  area 
S,  a  set  of  windows  {W^  ,W2,. . .Wn)  and  a  set  of  search  areas  (s\s2,..., 
Sn>  are  created  as  shown  in  Figure  2-4.  Resolution  of  the  window  and 
search  area  in  the  X  or  Y  direction  at  any  level  is  twice  the  resolution 
of  the  window  and  search  area  for  the  next  level,  respectively.  W1  can 
be  created  from  Wi-^  by  dividing  vf  into  blocks  of  size  2x2  and  treat¬ 
ing  each  block  as  a  pixel  with  value  equal  to  the  average  of  its  four 
pixels.  Similarly,  S1  can  be  created  from  S^.  Matching  starts  at  the 

lowest  resolution  level  (level -n)  where  Wn  of  size  —  x  —  is  matched 

~n  „n 

n  MM  t  c 

with  S  of  size  x  .  The  amount  of  computation  at  this  level, 

Z  2n 

therefore,  is  proportional  to  -  —■  +  1)(^-  -  ~  +  1).  Based 

2n  2n  2n  2n  2n  2n 

on  some  predetermined  criterion,  only  the  most  promising  test  locations 
are  selected  for  testing  in  (n-l)th  level. 


Figure  2-4.  Search  area  and  window  for  various  levels. 


In  the  (n-l)^  level,  Wn-^  is  matched  with  Sn”^  only  at  locations 
selected  in  nth  level.  This  procedure  continues  from  the  (n-1 level 
to  the  (n-2)tfl  level  and  so  on.  Finally  at  the  0th  level,  the  registra¬ 
tion  point  is  identified.  Savings  in  computation  depends  on  the  number 
of  levels  used  and  the  criteria  used  to  select  promising  test  locations 
at  each  level.  It  should  be  pointed  out,  however,  that  the  more  points 
eliminated  at  each  level,  the  greater  is  the  possibility  of  obtaining  a 
false  match. 

In  this  chapter  various  existing  methods  of  accomplishing  digital 
image  registration  are  presented.  The  most  commonly  used  method  is  cross¬ 
correlation.  There  exist  two  independent  ways  of  computing  the  correla¬ 
tion  surface  (Direct  method  and  Fast  fourier  transform  method).  The  FFT 
method  requires  a  large  amount  of  memory  for  software  implementation  and 
is  very  complex  for  real  time  hardware  implementation.  The  direct  method 
requires  less  memory  as  compared  to  the  FFT  method,  involves  no  complex 
multiplications  or  complex  additions,  and  can  be  easily  implemented  in 
real  time  using  digital  hardware.  Even  though  the  vector  correlation 
algorithm  is  expected  to  yield  better  performance  than  the  standard  cor¬ 
relation  algorithm,  its  use  is  limited  by  the  large  amount  of  computa¬ 
tion  required  to  implement  the  method  (more  than  twice  the  computation 
needed  by  the  standard  correlation  algorithm).  Due  to  the  difficulty 
encountered  in  the  automatic  segmentation  of  digital  images  into  homoge¬ 
neous  regions,  the  two  variations  of  the  standard  correlation  algorithm, 
namely,  feature  matching  correlation  and  hybrid  correlation  algorithms 
may  not  be  of  any  significant  use.  The  other  template  matching  method, 
sequential  similarity  detection  algorithm,  computes  an  error  surface  as 
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a  measure  of  dissimilarity  between  the  window  and  the  search  area. 


Since  this  method  requires  addition  and  a  few  division  operations,  it 
can  be  easily  implemented  (addition  and  subtraction  operations  are  sim¬ 
pler  than  multiplication  and  division  operations). 

Algorithms  which  accomplish  image  registration  by  matching  moments 
(moments  method)  or  straight  line  edge  content  of  the  window  and  the 
search  area  are  called  feature  matching  algorithms.  The  moments  method, 
which  is  computationally  inefficient  for  single  image  registration,  looks 
promising  for  multiple  image  registration.  The  necessity  of  a  highly 
efficient  algorithm  to  transform  a  digital  image  to  a  binary  image  (edge 
and  non-edge  pixels)  and  a  connectivity  test  to  identify  true  straight 
line  edges  (i.e.,  composed  of  adjacent  edge  pixels)  makes  the  use  of 
Hough's  transformation  for  digital  image  registration  less  attractive. 
Therefore,  it  is  concluded  that  the  standard  correlation  algorithm,  the 
SSDA  and  the  moments  method  are  more  promising  using  present  state-of- 
the-art  hardware.  A  detailed  comparison  of  the  above  three  methods  is 
presented  in  the  next  chapter.  It  should  be  pointed  out,  however,  that 
very  large  scale  integrated  circuits  developed  to  perform  a  specialized 
task  may  at  some  future  date  make  any  of  the  algorithms  discussed  above 
as  being  computationally  inefficient  feasible.  Because  of  this  possibil¬ 
ity,  a  follow-on  program  should  investigate  the  computational  accuracy 
of  some  of  these  methods. 


i  , 


III.  COMPARISON  OF  METHODS  FOR  MULTIPLE  IMAGE 
REGISTRATION 


Problem  of  Multiple  Image  Registration 
Let  the  search  area  S  and  n  windows  W-j ,  Wg,  . Wn  be  defined  as 
shown  in  Figure  3-1.  S  is  a  MxN  array  of  digital  picture  elements  which 
may  assume  one  of  G  possible  levels  on  the  gray  scale. 

0<S(1,j)<G-1  (3-1) 

for  1  <  j  <  N 

Wk  is  a  KxL  array  of  pixels  having  the  same  gray  scale  range. 

0  _<  VJkU,m)  <.  G-l  (3-2) 

for  1  <_  i  <_  K,  1  <  in  <  L  and  k  =  1 ,  2,  ....  n 

Each  KxL  subimage  of  S  can  be  uniquely  identified  by  its  upper  left  cor¬ 
ner's  coordinates.  Let  S.  .  denote  the  KxL  subimage  of  S  whose  upper 

*  *  J 

left  corner  is  (i,j). 

S,  ,(z,m)  =  S(i+t-l ,  j+m-1)  (3-3) 

1  >  J 

for  1  <  i  <  K,  1  <_  m  <_  L 

and  1  <  i  <  M-K+l,  1  _<  j  <  N-L+l 

If  S  and  do  not  differ  in  pixel  resolution  and  rotation,  the  multiple 

image  registration  problem  reduces  to  that  of  finding  (i£,  j£)  such  that 

Sj*  .*  best  matches  W.  ,  for  k  =  1 ,  2,  . . . ,  n. 
k’Jk 

Correlation  and  sequential  similarity  detection  algorithms  are 
commonly  used  for  the  registration  of  a  smaller  image  within  a  larger 
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image.  These  methods  have  been  proven  to  be  reliable  and  computation¬ 
ally  efficient  for  the  problem  of  single  image  registration.  However, 
a  method  which  is  efficient  for  single  image  registration  may  not  be 
efficient  for  multiple  image  registration.  In  general,  in  order  to  com¬ 
pare  the  computational  efficiency  of  different  algorithms,  the  number  of 
arithmetic  operations,  memory  requirement,  computational  speed  and  com¬ 
plexity  of  implementation  must  be  considered.  The  method  which  requires 
fewer  arithmetic  operations  for  software  implementation  may  be  complex 
or  even  infeasible  for  hardware  implementation.  It  is  extremely  difficult 
to  arrive  at  a  valid  means  of  comparison  without  knowing,  exactly,  how  the 
methods  are  implemented.  Therefore,  the  problem  of  multiple  image  regis¬ 
tration  is  studied  from  both  software  as  well  as  hardware  points  of  view, 
independently. 


Comparison  of  Software  Implementations 

If  the  algorithms  are  implemented  entirely  using  software,  the 
amount  of  core  memory  required  and  the  number  of  arithmetic  operations 
to  be  performed  can  be  used  as  a  means  of  comparison.  It  is  assumed 
that: 

a.  Multiplications,  divisions  and  squaring  operations  are 
equivalent  (i.e.,  1  multiplication  =  1  division  =  1 
squaring). 

b.  Enough  memory  is  available  to  store  the  intermediate 
results  for  future  use. 


c.  Computation  time  is  directly  proportional  to  the  number 
of  arithmetic  operations  performed  to  implement  the  method. 


Multiplication  and  addition  requirements  for  the  correlation,  the  SSDA 
and  the  moments  methods  are  derived  in  the  following  sections. 


Correlation  Algorithm 

Elements  of  normalized  correlation  function  of  search  area  S  and 
window  are  defined  to  be 


K  L  j 

t  X  X  WkU,m)  S.^U.m)] 

Rk(U)  =  x-i1 - g — [ - 


[I  I  W‘U,m)][X  X  sf  At. m)] 
£*1  m=l  *  1=1  m=l 


1  <  i  <  M-K+l  ,  1  1  j  l  N-L+l 
k  =  1 ,  2,  ...»  n 


.  K  L 

A;  :  s  [  I  X  w.(i,m)  S.  Ai,m)V 

1,J  i=l  m=l  *  1,J 

K  L  9 

h  *  Z  t  «£(*»■) 

*  m=l  K 


(3-4) 


(3-5) 


(3-6) 


K  L 


^■j  -j =  x  a  j(i»m) 
1,J  i=l  m=l  1,J 


(3-7) 


Therefore, 


A^  • 

^  _  _  .i>j _ 

k 

1  <  i  <  M-K+l  ,  1  <  j  <  N-L+l 


(3-8) 


k  *  1 ,  2,  . . . ,  n 
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Equations  for  the  number  of  multiplications  and  additions  in  terms  of 

M,  N,  K,  L  and  n  are  derived  in  the  four  steps  below: 

2  k 

1.  Computation  of  R.  (i,j)  from  A.  .,  8.  and  C.  .  requires 

K  \  9  J  K  I 

one  multiplication  and  one  division  or  two  equivalent 
multiplications.  Since  there  are  (M-K+l )( N-L+l )  ele¬ 
ments  in  each  of  the  n  correlation  surfaces,  a  total 
of  2n(M-K+l ){N-L+1 )  equivalent  multiplications  are  re¬ 
quired  for  this  stage. 

L 

2.  To  compute  A?  .,  (KL+1)  multiplications  and  (KL-1)  ad- 

'  *J 

ditions  are  performed.  Since,  for  each  of  the  n  win- 

i/ 

dows,  A.  .  must  be  computed  at  (M-K+l )(N-L+1 )  reference 

•  9  J 

points,  this  task  requires  a  total  of  (M-K+l )(N-L+1 ) (KL+1 )n 
multiplications  and  (M-K+l )(N-L+1 ) (KL-1 )n  additions. 

3.  Bk  is  computed  for  each  of  the  n  windows  and  thus  re¬ 
quires  KLn  multiplications  and  (KL-l)n  additions. 

4.  Cj  j  for  1*1,2,...,  M-K+l  and  j  =  1,  2,  ....  N-L+l 
can  be  computed  in  many  ways.  In  this  report,  it  is 
assumed  that  the  $T  .(£.,m)  values,  for  1  <  t  <  K 

'  9  J 

and  1  <  in  <  L,  are  first  computed  and  stored. 

Then  C.  *  is  computed  as 

*  *  J 


for  1  <  i  <  M-K+l  and  1  <_  j  <_  N-L+l 

Computation  of  C.  . 's  is  done  only  once  and  the  values  are  stored  in 

1  tJ 

memory  for  later  use.  The  above  tasks  require  MN  multiplications  and 
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(M-K+1 )(N-L+1 ) {KL-1 )  additions.  This  technique  requires  an  additional 
MN+(M-K+1 ) (N-L+l )  memory  locations. 

By  adding  the  partial  results  of  the  above  four  steps. 

Total  number  of  additions  =  (3-10) 

[(M-K+l) (N-L+l)  +  l][(KL-l)n]  +  (M-K+l ) (N-L+l ) (KL-1 ) 

Total  number  of  multiplications  *  (3-11) 

[ (M-K+l) (N-L+l )(KL+3)+KL]n  +  MN 

For  M  =  240,  N  =  256  ,  K  =  L  =  32, 

Total  number  of  additions  =  (3-12) 

48107598n  +  48106575 

Total  number  of  multiplications  *  (3-13) 

48295699n  +  61440 

Sequential  Similarity  Detection 
Algorithm' 

The  correlation  method  yields  the  correlation  surface  as  a  measure 
of  similarity,  while  the  sequential  similarity  method  computes  the  error 
surface  as  a  measure  of  dissimilarity.  The  normalized  error,  ek(i,j,£,m), 
between  W.(£,m)  and  S.  .(£,m)  is  defined  as 

K  1  t  J 

ek(i,j,£,m)  =  |SjjU»m)  -  -  WkU,m)  +  Wj  (3-14) 

where 

-  -I  K  L 

Wk  -  4-  1  I  W.(t,m)  (3-15) 

K  £■ 1  m=l  K 


and 
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1  K 
'*  A 


L 

E  Sj  ■  ( a , m } 
m=1 


(3-16) 


The  normalized  error,  E^C i , j ) ,  associated  with  the  reference  point  (i,j) 
is  defined  as 


K  L 

Ek(i , j )  *  S  Z  ek(i,j,i,m)  (3-17) 

ji=1  m=l 

for  1  <  i  <  M-K+l  ,  1  <.  j  <_  N-L+l 

To  register  Wk  within  S,  (i£,j£)  must  be  found  such  that 

Ek(ik*Jk)  <  for*  1  1  1  1  M-K+l.  1  t  1J  (3-18) 

and  1  <  j  <  N-L+l ,  j  t  j* 

Equations  for  the  number  of  multiplications  and  additions  required  to 
register  n  windows  using  the  SSOA  method  are  derived  in  the  following 
steps: 

1.  To  compute  each  Wk,  (KL-1)  additions  and  one  multipli¬ 
cation  (division)  are  required.  Since  there  are  n  such 
windows,  a  total  of  (KL-l)n  additions  and  n  multipli¬ 
cations  are  performed  to  compute  for  k  *  1 ,  2,  3, 

•  •  *  i  n  • 

A 

2.  It  is  assumed  that  S.  .  is  computed  independently  at 
each  reference  point  and  stored  in  memory  for  later 
use.  This  task  needs  (M-K+l ) (N-L+l ) (KL-1 )  additions 
and  (M-K+l ) (N-L+l )  multiplications  and  requires 
(M-K+l ) (N-L+l )  memory  locations. 

3.  Computation  of  Ek(1,j)  requires  (3KL)+(KL-1)  =  (4KL-1 ) 
additions.  Since  there  are  (M-K+l ) (N-L+l )  elements 


48 

In  each  of  the  n  error  surfaces,  a  total  of  (M-K+l) 

(N-L+l )(4KL-1 )n  additions  are  performed  to  compute 
n  error  surfaces. 

Therefore, 

Total  number  of  additions  *  (3-19) 

[(M-K+l )(N-L+1)(4KL-l)+(KL-1)]n  +  (M-K+l ) (N-L+l ) (KL-1 ) 

Total  number  of  multiplications  *  (3-20) 

(M-K+l) (N-L+l)  +  n 

For  M  *  240,  N  =  256,  K  =  L  *  32 

Total  number  of  additions  =  (3-21) 

1 92 5683 93 n  +  48106575 

Total  number  of  multiplications  =  (3-22) 

47025  +  n 

The  number  of  multiplications  required  is  negligible  when  compared  to 
the  number  of  additions. 

Moments  Method 

Given  a  two  dimensional  function  f(x,y),  the  moment  m  of  order 

pcj 

(p+q)  is  defined  by  the  relation, 

00  00 

V  ■  1  /  xpyp  f(x,y)dx  dy  (3-32) 

for  p,q  =  0,  1,  2,  3,  ... 

A  uniqueness  theorem  states  that  if  f(x,y)  is  piecewise  continuous  and  has 
non-zero  value  only  in  a  finite  region  of  X-Y  plane,  then  the  moments  of 
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all  order  exist  and  the  moment  sequence,  {11^},  1 s  uniquely  determined 
by  f(x,y)  and  conversely,  {n^}  uniquely  determines  f(x,y)  [10].  For 
a  digital  image  these  moments  are  given  by 


<"pq  x  1 I  xPyq  f(*.y) 

x  y 

for  p,q  *  0,  1,  2,  3,  ... 


(3-24) 


Since  the  digital  image  satisfies  all  the  conditions  required  by  the 
uniqueness  theorem  as  stated  above,  moments  of  all  order  exist  and  the 
moment  sequence  {mpcj>  can  be  used  to  accomplish  digital  image  registra 
tion  as  described  in  the  three  steps  below. 

1.  Let  {m  }  denote  the  moment  sequence  for  the  window 

pq 

W^.  Compute  all  moments  of  order  less  than  or  equal 
to  r  where  r  is  a  predetermined  number,  for  each  of 
the  n  windows. 


m*  »  I  I  Wk(i,m)  (3-25) 

fc=l  m=l 

for  p,q  *  0,  1 ,  2,  3,  . . . ,  r,  such  that  p+q  <_  r 

and  k  *  1 ,  2,  ...»  n. 

Although  the  reliability  of  this  method  increases  with 
increasing  r,  the  amount  of  computation  increases  with 
r  and  therefore  a  trade-off  exists.  For  some  applica¬ 
tions,  the  moments  method  has  been  found  to  be  success¬ 
ful  for  r  equal  two  or  three  [10]  -  [14]. 


2.  Compute  all  moments  of  order  r  and  less  for  each  of 
the  (M-K+l ) (N-L+l )  subimages  of  S. 
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for 

and 


p,q  =  0,  1,  2,  ....  r,  such  that  p+q  <_  r 
1  <  i  <  M-K+1,  1  <  j  <  N-L+l. 


(3-26) 


3.  To  register  Wk  within  S,  find  (i£,j£)  such  that 


for 

and 


m  (mpq 


-  <  l  (mk  -  m'J)2 

“/l  •  r\/-i  '  nn  • 


pq  pq  pq 

1  <  i  <  M-K+1  ,  i  +  1* 


pq' 


1  <,  j  <  N-L+l  ,  j  f 


(3-27) 


To  register  n  windows  W-j ,  Wg,  ....  Wp  within  S,  moments  for  the  n 
windows  and  (M-K+1 ) (N-L+l )  subimages  are  computed  only  once  (step  1  and 
2)  and  step  3  is  repeated  n  times  (once  for  each  window).  The  amount  of 

computation  associated  with  step  3  is  negligible  when  compared  to  the 

amount  of  computation  associated  with  steps  1  and  2.  Therefore,  although 
this  method  is  computationally  inefficient  for  n=l ,  its  efficiency  with 
respect  to  other  methods  increases  for  large  n. 

The  total  number  of  multiplications  and  additions  required  to  im¬ 
plement  this  method  depends  an  the  number  of  moments  used  to  characterize 
each  KxL  subimage.  Equations  for  the  number  of  additions  and  multipli¬ 
cations  required  to  register  n  windows  are  derived  for  the  following  two 
cases: 

a.  Case  1 :  All  moments  of  order  two  and  less  are  used  to  accomplish 

image  matching.  To  compute  these  six  moments  for  each  KxL  subimage, 

6(KL-1)  additions  and  8KL  multiplications  must  be  performed  (see  Table 
3-1).  There  are  n  windows  and  (M-K+1 ) (N-L+l )  subimages  of  size  KxL. 
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Table  3-1 


Computation  of  moments  of  order 

three  and  less. 

Moment 

Additions 

Multi  pi ications 

o 

O 

B 

(KL-1 ) (M-K+l ) (N-L+l ) 

0 

m01 

(KL-1 )(M-K+1 )(N-L+1 ) 

KL (M-K+l) (N-L+l) 

m10 

(KL-1) (M-K+l) (N-L+l) 

KL  (M-K+l)  (N-L+l) 

mn 

(KL-1) (M-K+l) (N-L+l) 

2KL (M-K+l) (N-L+l) 

m02 

(KL-1) (M-K+l) (N-L+l) 

2KL(M-K+1) (N-L+l) 

m20 

(KL->) (M-K+l) (N-L+l) 

2KL (M-K+l) (N-L+l) 

m12 

(KL-1) (M-K+l) (N-L+l) 

3KL (M-K+l) (N-L+l) 

m21 

(KL-1) (M-K+l) (N-L+l) 

3KL(M-K+1) (N-L+l) 

m30 

(KL-1) (M-K+l) (N-L+l) 

3KL (M-K+l) (N-L+l) 

m03 

(KL-1) (M-K+l) (N-L+l) 

3KL(M-K+1) (N-L+l) 
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To  register  each  of  n  windows  within  S,  11 (M-K+l )(N-L+1 )  additions  and 
6(M-K+1 )(N-L+1 )  multiplications  must  be  performed  (step  3).  Therefore, 

Total  number  of  additions  =  (3-28) 

6(KL-1) (M-K+l )(N-L+1)  +  [11 (M-K+l ) (N-L+l )  +  6(KL-l)]n 

Total  number  of  multiplications  =  (3-29) 

8KL(M-K+1) (N-L+l)  +  [8KL  +  6 (M-K+l ) (N-L+l )]n 

For  M  =  240,  N  =  256,  K  ■  32  and  L  =  32 

Total  number  of  additions  =  (3-30) 

288639450  +  52341 3n 

Total  number  of  multiplications  =  (3-31) 

385228800  +  290342n 

b.  Case  2:  All  moments  of  order  three  and  less  are  used  to  accomplish 
image  matching.  Following  the  steps  of  case  1  it  can  be  shown  that 

Total  number  of  additions  =  (3-32) 

481065750  +  903705n 

Total  number  of  multiplications  *  (3-33) 

963072000  +  490730n 

Comparison  of  Multipl ication  and 
Addition  Requirements 

The  number  of  additions  and  multiplications  which  are  required  to 
implement  correlation,  SSDA  and  moments  method  for  various  values  of  n, 
are  shown  in  Figures  3-2  and  3-3,  respectively.  From  Figures  3-2  and 
3-3  it  is  clear  that  if  n  is  greater  than  7,  moments  method  with  r  equal 
two  requires  less  computation  when  compared  to  the  correlation  method. 


additions  to  register  n  windows 


Figure  3-3.  Number  of  multiplications  to  register  n  windows 
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However,  if  r  equals  3,  n  must  be  greater  than  19  for  the  moments  method 
to  have  a  computational  advantage  over  the  correlation  method.  Substan¬ 
tial  savings  in  computation  can  be  accomplished  for  the  moments  method 
by  having  lookup  tables  for  each  of  the  moments  as  described  below. 

The  moment  mpq  of  a  two-dimensional  discrete  function  f(x,y)  of  size 
KxL  is  given  by 

K  L  n  n 

m  *  Z  J  *Pmq  fU,m)  (3-34) 

2=1  m=l 

2,Pmq  is  a  constant  for  given  values  of  i,  m,  p  and  q. 

Let,  Km  =  ^  (3-35) 

Now,  if  for  t  *  1,  2,  3,  ...,  K  and  m  =  1,  2 . L  are  precomputed 

and  stored  in  memory,  mpq  can  be  computed  by  performing,  only,  KL  multi¬ 
plications  and  (KL-1)  additions  no  matter  what  the  values  of  p  and  q  are. 
Assuming  that  such  look-up  tables  are  available  for  all  moments  of  order 
3  and  less,  equations  for  the  number  of  additions  and  multiplications 
are  recomputed. 

For  r  equal  two; 

Number  of  additions  =  (3-36) 

6 ( KL-1 ) (M-K+l ) (N-l+1 )+[l 1 (M-K+l ) (N-L+l )  +  6(KL-1 )]n 

Number  of  multiplications  = 

5KL(M-K+1 )(N-L+1#)+[5KL+6(M-K+1 ) (N-L+l )]n  (3-37) 

By  substituting  values  for  M,  N,  K  and  L, 

Number  of  additions  =  (3-38) 


288639450  +  52341 3n 
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Number  of  multiplications  *  (3-39) 

240768000  +  287270n 

Similarly,  if  r  equal  three. 

Number  of  additions  =  (3-40) 

481065750  +  9037 05n 

Number  of  multiplications  =  (3-41) 

433382400  +  479466n 

The  number  of  additions  remained  the  same,  but  the  number  of  multi¬ 
plications  is  reduced  substantially.  This  fact  is  graphically  shown,  in 
dotted  lines,  in  Figure  3-3.  In  general,  an  accurate  comparison  of  com¬ 
putation  time  required  by  different  methods  is  not  possible.  This  is 
because  the  ratio  of  multiplication  time  to  addition  time  depends  on  a 
number  of  factors  such  as  machine  used,  bit  length  and  algorithm  used  to 
accomplish  multiplication  of  two  numbers.  However,  the  previous  analysis 
shows  that  the  moments  method  is  computationally  feasible  and  takes  less 
number  of  arithmetic  operations  than  the  correlation  method  if  the  number 
of  windows  to  be  registered  is  sufficiently  large.  A  direct  comparison 
of  the  SSDA  and  moments  methods  is  difficult  because  the  SSDA  method 
requires  more  additions  while  the  moments  method  requires  more  multipli¬ 
cations.  Since  the  ratio  of  the  time  required  to  perform  one  real  multi¬ 
plication  to  the  time  required  to  perform  one  real  addition  is  not  known, 
one  single  measure  of  comparison  cannot  be  determined.  In  general, 
real  multiplication  time  *  a  x  real  addition  time, 
where  the  value  of  "a"  depends  on  the  machine,  its  bit  length  and  algo¬ 
rithm  used  to  implement  multiplication.  It  is  bounded  by  the  bit  length. 
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however  (i.e. ,  a  <  bit  length).  Assuming  that  'a'  equals  three,  for  a 
machine  such  as  the  IBM  370  the  total  number  of  equivalent  real  addi¬ 
tions  required  to  implement  correlation,  SSDA  and  moments  methods  (for 
r  equal  2  and  3)  for  different  values  of  n  are  plotted  in  Figure  3-4. 
If  all  moments  of  order  two  or  less  are  used  to  characterize  each  KxL 
subarray,  the  moments  method  using  lookup  tables  is  computationally 
more  efficient  than  the  other  methods,  for  n  greater  than  or  equal  to 
five. 


Comparison  for  Hardware  Implementation 
To  accomplish  multiple  image  registration  in  real  time  (or  almost 
so),  image  matching  algorithms  must  be  implemented  using  fast  hardware. 

A  given  algorithm  can  be  implemented  in  many  ways.  Three  simple  sche¬ 
matics  shown  in  Figure  3-5,  3-6,  3-7,  one  for  each  method,  are  used  for 
comparison  of  the  complexity  of  implementation. 

Correlation  Algorithm 

Correlation  between  the  window  W.  and  the  subimage  S-  .  is  com- 

K  •  *J 

puted  in  four  stages  as  described  below. 

2  2 

Stage  1:  In  the  first  stage,  S.  .(n,m),  Wr(4,m)  and  the  product 

1  >  J  K 

S.  . U,m)W.  (£,m)  are  computed  for  i  -  1,  2,  3,  . ...  K  and  m  =  1, 

I  >  J  K 

2,  ...,  L.  This  requires  a  total  of  3KL  two-input  multipliers 
(three  multipliers  for  each  pixel  pair). 

rr  K  l 

Stage  2:  In  the  second  stage,  If  A?  .  *  T  Y  S.  *  U,m)W,  (i,m), 

4=1  m=l  1,J  k 

K  L  «  K  l 

=  Z  Z  wl. ( )  and  C.  ,  =  Y,  Z  .-(t,m)  are  computed 

K  t-l  m-1  k  1,J  a=l  m=l  1>J 


using  three  KL-input  adders. 


Multiplier 


3-5.  Schematic  for  correlation  method. 


3-6.  Schematic  for  sequential  similarity  detection  algorithm. 


(i,,m)  S.  <(K»L) 


Figure  3-7.  Schematic  for  the  moments  method. 
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Stage  3:  In  the  third  stage,  two  two-input  multipliers  compute 

U 

A?  .  and  the  product  B.C.  .. 

1  h  IjJ 

2 

Stage  4:  Finally,  a  divider  computes  R£(i,j)  using  the  results 
obtained  in  Stage  3. 

Therefore,  3KL+3  two-input  multipliers  and  three  KL-input  adders  are  re- 

2 

quired  to  compute  Rf(i,j)  from  S.  .  and  W. . 

K  I )  J  K 

Sequential  Similarity  Detection 
Algorithms 

The  error  between  the  window  W,  and  subimage  S.  .is  computed  in 

K  1  ,  J 

three  stages  as  described  below. 

Stage  1:  In  the  first  stage,  $.  .  and  W.  are  computed  using  two 

1  »  J  K 

KL-input  adders  and  two  dividers. 

Stage  2:  In  the  second  stage,  ek(i,j,£,m)  for  i  *  1,  2,  ....  K 
and  m  s  1,  2,  ...»  L  are  computed  using  KL  four-input  adders. 

Stage  3:  Finally,  one  KL-input  adder  computes  Ek(i,j)  as  shown 
in  Figure  3-6. 

Therefore,  SSDA  implementation  requires  three  KL-input  adders,  KL  four- 
input  adders  and  two  multipliers. 

Moments  Method 

A  schematic  for  the  computation  of  all  moments  of  order  two  and 

less  is  shown  in  Figure  3-7.  Moments  of  S,  .  can  be  computed  in  two 

* 

stages. 

Stage  1:  In  the  first  stage,  the  products  IS.  -U.m),  mS.  .(t,m), 

2  2  1,J 
unS.  i  S.  .(i,m)  and  nrS.  .(s,,m)  are  computed  for  i  -  1 , 

*  »  J  *  I  »  J 

2,  ....  K  and  m  =  1,  ....  L.  5KL  two-input  multipliers  are  re¬ 


quired  for  the  above  purpose. 
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Stage  2:  In  the  second  stage,  six  KL-input  adders  compute  rrtQQ, 

m10*  m01’  mll’  m20*  and  m02‘ 

Therefore,  the  moments  method  requires  a  total  of  4KL-5(K+6)+6  two- 
input  multipliers  and  six  KL-input  adders  to  compute  all  moments  of 
order  two  and  less. 

From  Figures  3-5,  3-6,  and  3-7,  it  can  be  seen  that  all  three 
methods  can  be  easily  implemented  using  multipliers  and  adders  (i.e., 
none  of  the  methods  involve  function  evaluation).  The  moments  method 
looks  simpler  than  the  other  two  methods  with  just  one  level  of  multipli¬ 
cations  and  one  level  of  additions.  In  general,  no  method  has  any  sig¬ 
nificant  advantage  over  others  as  far  as  the  complexity  of  implementation 
is  concerned. 

If  the  video  images  are  sampled  at  5  MHz  and  if  it  is  required  to 

register  all  n  windows  simultaneously  (or  almost  so)  between  sampling, 

more  than  one  correlator  (probably  n  correlators)  may  have  to  be  used  in 

parallel  when  implementing  the  correlation  or  SSDA  method.  Therefore,  if 

n  is  large,  these  two  methods  may  prove  uneconomical.  The  maximum  number 

of  parallel  units  that  can  be  used  is  also  restricted  by  cost  and  the 

volume  of  space  and  power  available.  For  the  moments  method,  however,  if 

there  is  one  unit  of  hardware  to  compute  all  the  required  moments,  the 

degree  of  mismatch  between  S.  .  and  each  of  the  n  windows  can  be  computed 

almost  simultaneously  as  shown  in  Figure  3-8.  The  hardware  required  to 

compute  £  (m  Jj  -  m^)^  for  k  *  1,  2,  ...,  n  is  negligible  and  therefore 

pq  pq  pq 

moment's  method  accomplishes  multiple  image  registration  using  relatively 
less  hardware  as  compared  to  the  correlation  and  SSDA  methods  when  n  is 


Figure  3-8.  Parallel  computation  of  dissimilarity  between  S.  .  and  each  of  the  n  windows. 
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In  this  chapter,  computational  efficiency  of  the  correlation, 
the  SSDA  and  the  moments  method  is  compared  from  software  as  well  as 
hardware  points  of  view,  independently.  It  is  found  that  the  moments 
method  becomes  more  efficient  as  the  number  of  windows  increases.  It 
takes  less  computation  time  if  implemented  using  software  and  less  hard¬ 
ware  for  real  time  implementation  if  the  number  of  windows  is  sufficient¬ 
ly  large.  In  general,  feature  matching  algorithms  are  expected  to  be 
more  efficient  than  template  matching  algorithms  when  the  number  of  win¬ 
dows  is  large.  Modifications  of  the  moments  method  and  two  new  feature 
matching  algorithms  are  presented  in  the  next  chapter. 
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IV.  NEW  METHODS  FOR  DIGITAL  IMAGE  REGISTRATION 

For  multiple  image  registration,  it  has  been  shown  that  the  mo¬ 
ments  method  is  computationally  more  efficient  than  template  matching  al¬ 
gorithms  if  the  number  of  windows  is  sufficiently  large.  This  is  because 
computation  for  template  matching  algorithms  is  directly  proportional  to 
the  number  of  windows  whereas  in  feature  matching  algorithms  features  are 
extracted  for  all  subimages  of  the  search  area  and  windows  only  once  and 
the  matching  procedure  is  repeated  once  for  each  window.  Computation  re¬ 
quired  to  match  the  features  is  negligible  compared  to  that  required  to 
compute  the  features. 

In  general,  feature  matching  algorithms  for  multiple  image  regis¬ 
tration  are  expected  to  be  more  efficient  than  template  matching  algo¬ 
rithms.  For  this  reason,  an  effort  vas  made  to  improve  the  moments  meth¬ 
od  and  to  develop  new  feature  matching  algorithms.  Two  new  feature 
matching  algorithms,  one  based  on  intraset  and  interset  distances  [26], 
and  the  other  based  on  correlation  between  adjacent  pixels  [27],  are 
presented  in  this  chapter.  The  computational  efficiency  of  each  of  the 
new  methods  is  compared  with  that  of  existing  methods.  A  technique  of 
making  the  standard  correlation  algorithm  more  suitable  for  multiple  reg¬ 
istration  is  also  described. 

Moments  Method 

In  order  to  obtain  meaningful  results  from  any  of  the  image  reg¬ 
istration  algorithms,  it  is  necessary  to  preprocess  windows  and  subimages 
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of  the  search  area  such  that  the  mean  and  the  standard  deviation  of  their 
pixel  values  are  equal.  This  is  called  intensity  level  normalization  and 
is  required  when  the  window  and  the  search  area  are  obtained  from  sensors 
with  different  d.c.  gain  and  bias.  Ideally,  each  KxL  subarray  of  the 
search  image  should  be  preprocessed  to  have  zero  mean  and  unity  standard 
deviation.  This,  however,  would  require  too  much  additional  computation. 
Fortunately,  normalization  can  be  incorporated  into  the  moments  method 
with  almost  no  additional  computation.  Two  cases  where  it  is  not  possi¬ 
ble  to  use  the  moment  sequence  {nip^ >  directly  to  accomplish  image  regis¬ 
tration  without  intensity  level  normalization  are  given  below.  Means  of 
incorporating  normal ization  within  the  algorithm  without  actually  chang¬ 
ing  the  pixel  values  are  also  described. 

Case  1:  Let  W  and  $.*  .*  be  the  window  and  its  matching  subimage 

I  f  J 

of  the  search  area,  respectively.  Since  the  two  images  are  obtained  from 
different  sensors,  the  corresponding  pixels  may  not  have  the  same  pixel 
value.  Assume  that  W  and  S.*  •*  are  related  by  the  following  equation. 

*  9  J 


W(z,m)  =  c  S-*  .*(t,m)  ,l<t<K,  l<m<L 

*  9  J  ~~ 


(4-1) 


where  c  is  a  constant. 

u 

The  moment,  m  ,  of  the  window  is  given  by 


=  Z  E  *Pmq  W(*,m) 
pq  z=l  m=l 


(4-2) 


for  p,q  =  0,  1,  2,  ... 

S  j  *  -j  * 

The  moment,  rn  ,  of  the  subimage  $.+  .*  is  given  by 

pq 
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K 

I 


L 

l 


S. 


(4-3) 


From  Equations  (4-1),  (4-2)  and  (4-3),  it  is  clear  that 


W  i*  i* 

rn  =  c  m  ‘  for  p,q 

pq  pq 


0,  1,  2,  ... 


U 

Any  attempt  to  match  the  moment  sequences  (m^ 
leads  to  false  registration.  However, 


}  and  {m  }  directly 

pq 


W  ^i*  i* 
-  < 

”“o  *0j*>J* 


for  p,q  =  0,  1 ,  2,  ... 


(4-4) 


m“_  m  1 -J 


Therefore,  the  sequences  {-^9}  and  {-£9 — >  must  be  matched  to  accomplish 

m00  m  1  » ^ 

registration.  uu  m00 

Case  2:  Assume  that  S.*  .*  and  W  are  related  as  described  by 
Equation  (4-5). 


W(t,m)  =  c  S-*  -*U, m)  +  d  ,  1  <  i  <  K  and  1  <  m  <  L  (4-5) 

1  » J  —  — 

where  c  and  d  are  constants. 

For  this  case, 

S . 

mpq  t  for’  P»q  s  0.  1.  2,  ...  (4-6) 

Let  W  and  S-*  .*  be  the  mean  pixel  values  of  W  and  S.+  .*,  respectively. 

1  l  ,j 

Suppose,  W  is  obtained  from  W  by  subtracting  its  mean  pixel  value  W 
from  each  pixel.  S'.*  .*  is  obtained  from  S.*  .*,  similarly.  Now,  W  and 

1  »j  i  ij 

S!*  4*  are  related  by  the  following  equation. 

*  f  J 

W'U,m)  =  c  S'  .*  (i,m)  (4-7) 

1  5  J 

1  <  i  <  K  and  1  <  m  <  L 


for 
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The  following  relations  between  the  moments  of  W  and  S!*  iic  can  be  eas- 
ily  derived 


moo  =  moo 


=  0 


and 
Therefore, 


U*  •?* 

fflpq  =  c  mpq  *  ,  for  p,q  =  0,  1,  2,  . . 


(4-8) 

(4-9) 


W' 

PQ  =  pq 
VF  S' 

ml  „  i*,J* 


rs  m 


rs 


for 


(4-10) 


p,q  —  0,  1,  2,  .  • . 


provided  that  r  and  s  are  not  simultaneously  equal  to  zero.  Therefore, 


W 


S’.  • 


m”  m 

for  this  case  the  sequences  l-fp4  and  {-§? — }  can  be  matched  to  accom- 

z  & 


plish  image  registration.  However, 


K  L  n 

W  Y  Z  *Pmq  (W(a,m)  -  W) 

mpq  _  £*1  m=l 

y*  K  l 

mrs  £  Z  s-rmS  (W(a,m)  -  W) 

i=  1  m=l 

Y  Y  *Pmq  W(t,m)  -  W  Y  Z  lPm<1 

_  £=1  m=l  _ a=l  m=l 


K  L 


K  L 


£  £  arms  W(t,m)  -  W  £  £ 

:=1  m=l  a=l  m=l 


„r  s 
a  m 


(4-11) 


> 
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K  L 


m 

=  _ea 


.  u  I  2  1 

pq  £.=1  m=l 


Pmq 


u  -  K  L  „  e 

m“  -  «  I  I  iV 

rs  £=1  m-1 


Now, 


.  i  K  L 

w  =  I  I  wu,m) 

KL  t-l  m=1 


_w 

m00 

KL 


(4-12) 


Let  K  and  K  be  constants  defined  as 
pq  rs 

K  =  i_  y  V  *Pmq 
Pq  KL  A  m=l 


(4-13) 


K  L 


and 


1  ^  _ 


(4-14) 


From  Equations  (4-11)  through  (4-14) 
mW'  mW  *  mW 

Pq  _  mpq  Kpqm00 

jr~  Jr  fjr 

mrs  mrs  '  Krsm00 


(4-15) 


K  and  K  „  are  constants  which  can  be  precomputed  and  stored.  Therefore, 
pq  rs 

w  „  mw 

m  —  K  jTOqq 

it  is  clear  from  Equation  (4-15)  that  if  sequences  {- -jp - )  and 

mrs  '  Krsm00 

Si  J  Si  J 

mDO  J  ■  Knam0Q  J 

(_2g - £2_uu — j  are  useC|  for  image  registration,  intensity  level 


1 »J  _  /  m  1 ’3 
rs  rsm00 


normalization  is  accomplished  with  almost  no  additional  computation. 

This  makes  the  moments  method  more  reliable  without  excessive  additional 
computation. 
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Distance  Measures  for  Digital 
Image  Registration 

Pattern  classification  by  distance  measures  is  one  of  the  earliest 
concepts  in  automatic  pattern  recognition.  The  motivation  for  using  dis¬ 
tance  functions  as  a  classification  tool  follows  from  the  fact  that  the 
most  obvious  way  of  establishing  a  measure  of  similarity  between  pattern 
vectors,  which  can  also  be  considered  as  points  in  Euclidean  space,  is 
by  determining  their  proximity.  Common  distance  measures  used  for  pre¬ 
processing  and  feature  extraction,  a  new  registration  method  based  on 
intraset  and  interset  distances,  the  number  of  arithmetic  operations  re¬ 
quired  to  implement  this  method  and  the  relation  between  distance  measures 
and  moments  of  a  digital  image  are  described  in  this  section. 

Distance  Measures 


Point-to-Point  Distance.  Let  P.  and  P.  be  any  two  points  in  the 
two-dimensional  x-y  plane.  Each  point  can  be  uniquely  identified  by  its 
x  and  y  co-ordinates,  i.e., 


pj 


xi 

yj 


The  square  of  the  distance  between  points  P.  and  P.  is  given  by 

J 

D2(prPj)  =  (Vxj)2  +  (YYj)2 


(4-16) 
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Point-to-set-di stance.  Let  P_  be  a  point  and  S,  be  a  set  of  'a* 
■  -  — ■  1 '  ~  -  -  —  . —  0  » 

points  in  the  two-dimensional  plane,  i.e.. 


The  distance  between  the  point  P  and  set  S_,  D^(P_ ,S  ) »  is  defined  as 

0  a  U  a 

the  mean  square  distance  between  the  point  PQ  and  the  *a'  members  of  the 
set  Sa.  Therefore, 

Q 


°2<P„>V  *  r  I,  [<VV2  +  <W2J 


(4-17) 


Intraset  distance.  The  intraset  distance  of  a  set  is  defined  as 
the  mean  square  distance  between  the  points  of  that  set  [26],  The  mean 
square  distance,  02(Pj ,Sfl),  between  a  fixed  point  P^  and  all  other  a-1 
points  of  the  set  Sa  is  given  by 

u 

D2(pj’sa}  =  irr  i  [(xrxi)2  +  (yV2J  (4_18) 

Since  D2(P,,P.)  is  zero,  D2(P.,Sj  can  be  written  as 
J  J  j  « 

°2<Pj'Sa'  *  HT  C(XrX()Za  (Yj-Y,)2]  (4-19) 

2 

Therefore,  the  intraset  distance  of  S  .  D  (Sj,  is  given  by 

a  a 

°2<sa>  -  5T7TT  j|,  C(xj-Xi)2  +  <W2’ 


(4-20) 
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Equation  (4-20)  can  be  easily  reduced  to  a  simple  closed  form  in  terms 
of  statistical  properties  as  described  below: 


02(Sa) 


•  i^TT  J,  I  ( YV2  'J&i  1 1, 


a  a 


(4-21 ) 


In  Equation  (4-21)  (Xj-X,)2  and  ^  £  (Yj-Y, 


a  a 


r 


will  be  referred  to  as  the  X-component  and  Y-component  of  intraset  dis¬ 
tances,  respectively. 

Let 


X  =  l  2  Xi  and  Y  -  £  £  Y  . 

a  i * 1  1  a  i*l  1 


w  _  1 


(4-22) 


Then 


a  a 


a  a 


_  _  2 

V  \4  l  V\fc 


I  I,  (Xj-X,)  -  I  X,  (Xj-X-X  +X) 
J=1  1=1  J  jil  1=1  J 


a  a 


=  I  l  [(xr*r  +  ^xi 

j=l  i=l  3  1 


2  .  /v  X)2 


2(Xj-I)(X.-X) 

(4-23) 


It  can  be  shown  that 


a  a 


T  I  (x.-x)(x.-x) 

j=l  i=l  J  1 


=  0 


(4-24) 


Therefore, 


3  3  a  3  3  «  3  3 

l  1  (Xrx.)2=  l  l  (X.-X)2+  l  I  (X.- 
j=l  i=l  J  1  j=l  isl  3  j=l  i=l  1 


X)2 


(4-25) 
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=  a 


£  tyx)2  *  a  £  (xr7): 


2  2  2  2 
=  a  ax  +  a  cx 


2  2 
"  2a  °x 


where 


<r  2  *  \  t  (X.-X)2  *  [  |  (X  -X)2  (4-26) 

x  a  i=l  1  a  j=l  J 

2 

Therefore,  the  X-component  of  intraset  distance,  Dv(S. ),  is  given  by 

X  a 

°x>sa)  *  rr  °x2  <4-27> 

Similarly,  the  Y-component  of  intraset  distance  is  given  by 

°y<Sa>  *  i?T  »y2  <4-28> 

where 

o  2  -  j  Z  (V,-T)2  (4-29) 

y  a  i=l  1 

From  equations  (4-2G),  (4-27)  and  (4-28), 

D2(Sa)  *  (ax2  +  Oy2)  (4-30) 


Equation  (4-30)  expresses  intraset  distance  in  terms  of  variances  associ¬ 
ated  with  the  X  and  Y  coordinates  of  the  set  points. 
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Interset  distance.  Let  Sa  be  a  set  of  'a'  points  and  Sb  be  a  set 

of  1 b '  points  in  the  two-dimensional  plane.  Then,  the  interset  distance 

between  sets  S  and  S.  is  defined  as  the  distance  between  their  centroids 
a  o 

[26].  If  (Xa«  7fl)  and  (Xb>  7b)  are  the  centroids  of  set  Sa  and  Sb>  re¬ 
spectively,  the  interset  distance  Squared  between  Sa  and  Sb  is  given  by 

°2<W  ■  <W2  +  <V7b>2  (4-3,) 


Feature  Vector  for  Digital  Image 

Let  F  be  a  KxL  array  of  digital  pixels  which  may  assume  one  of  G 


gray  levels. 


1  <  F(X,Y)  <G,1<)(<K,  1<Y<L  (4-32) 

Let  S9  denote  the  set  of  pixels  with  pixel  value  equal  to  g  and  ng  be  the 

number  of  elements  in  S9.  Now,  digital  image  F  can  be  considered  as  the 

12  C 

union  of  sets  S  ,  S  ,  ...,  S  . 


F  =  S1  U  S2  U  ...  U  SG  (4-33) 

Each  element  of  the  set  S9,  which  is  just  a  point  in  the  two-dimensional 
X-Y  plane,  is  uniquely  specified  by  its  coordinates. 


where 


(4-34) 


(4-35) 


and 


g  =  1,  2, 


•  •  •  » 


G 
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Each  of  the  6  sets  can  be  characterized  by  its  intraset  distance  which  is 
defined  as  the  mean  square  distance  between  its  elements.  From  Equation 
4-3  the  intraset  distance  of  the  set  S9  is  given  by 


where 


D2(S9)  =  D2(S9)  +  D2(S9) 


and 


1 


(Y?-?9)2 


(4-36) 


(4-37) 


(4-38) 


The  location  of  one  set  with  respect  to  the  other  can  be  characterized  by 
the  interset  distance  between  them  which  is  defined  as  the  distance  be¬ 
tween  their  centroids.  Since  there  are  G  sets  in  F,  there  exist 
distinct  interset  and  G  intraset  distances.  A  ^  -dimensional  feature 
vector  consisting  of  G  intraset  and  -G(G— )  interset  distances  maps  the 
digital  image  F  to  a  point  in  the  -dimensional  Euclidean  space. 


Derivation  of  Multiplication 
and  Addition  Requirements 

Consider  the  digital  image  F  of  size  KxL  described  by  Equations 
(4-32)  and  (4-33).  Intraset  distance  of  the  set  S9  is  given  by  Equation 
(4-36).  The  X-component  of  the  intraset  distance,  D2(S9),  can  be  written 


as 
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(4-39) 


where 


-  s9  A 

g  1=1  1 


(4-40) 


The  multiplication  and  addition  requirements  to  compute  the  feature 
vector  are  derived  in  the  steps  below.  It  is  assumed  that  multiplica¬ 
tions,  divisions  and  squaring  operations  are  equivalent. 

1.  From  Equation  (4-40),  it  is  clear  that  n^-1  additions  and 
one  multiplication  (division)  are  required  to  compute  P. 

2.  Computation  of  D^(S^)  as  described  by  Equation  (4-39)  re¬ 
quires  n^+4  multiplications  and  n^+1  additions.  Therefore, 

a  total  of  2n  additions  and  (n  +5)  multiplications  are 
9  9 

needed  for  the  computation  of  D^(S^). 

3.  Similarly,  2n  additions  and  (n  +5)  multiplications  are 
required  for  the  computation  of  D^(S^). 

4.  Adding  results  of  steps  2  and  3,  4nQ  additions  and  2nQ  +  10) 
multiplications  must  be  performed  to  compute  D^(S^). 
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5.  Since  there  are  G  sets  in  the  image  F, 

G 

Total  number  of  additions  =  £  (4n  )  =  4KL  (&-41 ) 

9=1  9 

G 

Total  number  of  multiplications  =  £  (2n  +10)=2KL+10G  (4-42) 

g=l  9 

6.  The  interset  distance  between  two  sets  is  the  distance 
between  their  centroids  (See  equation  (4-31)).  Two  multi¬ 
plications  and  three  additions  are  required  to  compute 

each  of  the  interset  distances. 

Total  number  of  additions  =  — (4-43) 

Total  number  of  multiplications  =  G(G-1 )  (4-44) 

Therefore,  in  order  to  compute  the  feature  vector  for  a  KxL  array  of  pix¬ 
els  which  can  assume  one  of  G  gray  levels,  a  total  of  4KL+  3GlG-~--) 
additions  and  2KL+G(G+9)  multiplications  must  be  performed. 

Multiple  Image  Registration 

A  method  of  accomplishing  multiple  image  registration  using  feature 

vectors  of  windows  and  subimages  of  S  is  described  below. 

Step  1:  Compute  the  feature  vector  for  each  of  the  n  windows  and  each  of 

(M-K+l ) (N-L+l )  subimages  of  S.  Let,  and  V1’3  denote  feature  vectors 

of  the  window  W,  and  subimage  S.  .,  respectively. 

K  I  9  J 
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Vk  = 


'1 

,k 


6(6+1 )  , 
2  J 


for  k  =  1 ,  2 ,  ....  n 


(4-45) 


.i.J 


for,  1  <_  i  <_  M-K+l 
and,  1  <  j  <  N-L+l 


(4-46) 


G(G+1) 

2  -> 

The  above  vectors  can  be  considered  as  points  in  -dimensional 

Euclidean  space. 

Step  2:  To  register  W.  within  S,  find  the  subimage  S.*  ■*  whose  feature 

VJk 

vector  best  matches  that  of  W.  .  In  other  words,  find  S.*  •*  such  that 

*  VJk 


i *  i* 

,k’Jk,  1 2 


^  _  yi * J  I  I  2 


(4-47) 


for  1  _<  i  _<  M-K+l,  i^i* 

and  1  js  j  <.  N-L+l , 

where,  |  [ V^-v"*  |  | ^  is  the  square  of  the  distance  between  vectors  and 

yi  ,3  ' 


.J 
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To  register  n  windows  W-j ,  . Wn  within  S,  the  feature  vec¬ 

tors  for  n  windows  and  (M-K+l ) (N-L+l )  subimages  are  computed  only  once 
(step  1)  and  step  2  is  repeated  n  times,  once  for  each  window.  The  com¬ 
putation  associated  with  step  2  is  small  when  compared  to  the  computation 
associated  with  step  1.  Therefore,  although  this  method  is  computation¬ 
ally  inefficient  for  n=l ,  its  efficiency  with  respect  to  the  correlation 
method  increases  for  large  n.  Equations  for  the  number  of  additions  and 
multiplications  required  to  register  n  windows  are  derived  below. 

In  step  1  (M-K+l ) (N-L+l }+n  feature  vectors  are  computed. 

4KL+3£l|lil 

additions  and  2KL+G(G+9)  multiplications  are  required  to  com¬ 
pute  one  feature  vector.  Therefore, 

Number  ofadd1t1ons}  =  [4KL  +  M|in][(M.K+1 )  {N.L+1  )+n]  (4-48) 

^cations^  i  rTstep^ 1  ~ *  =  [2KL  +  G(G+9)][(M-K+1 ) (N-L+l  )+n]  (4-49) 

In  step  2,  G2+G-l  additions  and  multiplications  are  performed  to 

compute  | | -  V^,J'||2.  There  are  (M-K+l ) (N-L+l )  reference  points  and  n 
windows.  Therefore, 


“10ns>  =  (G2+G-l)  (M-K+l)  (N-L+l  )n  (4-50) 

Mimber  of  multipli-,  _  G ( G+l )  ,  .,N  ,n  C1  ■ 

cations  in  step  2  1  2~^  (M-K+l ) (N-L+l )n  (4-51, 

From  Equations  (4-48)  through  (4-51), 


Total  number  of  additions  = 


-K+l) (N-L+l )+n] 


+ 


(G2+G-l  )(M-K+1 ) (N-L+l )n 


(4-52) 


(4-53) 


82 

Total  number  of  multiplications  = 

[2KL+G(G+9)][(M-K+1 )(N-L+1 )+n]  +  G(G±ll  (M-K+l ) (N-L+l )n 
For  M  =  240,  N  =  256  and  K  =  L  =  32, 

Total  number  of  additions  =  (4-54) 

[47025  (G2-t-G-l)  +  1.5  G(G-l)  +  4096]n  +  70537.5  G(G-1 )  + 

192614400 

Total  number  of  multiplications  =  (4-55) 

[2048  +  G(G+9)  +  23512. 5(G+l)G]n  +  47025  G(G+9)  +  96307200 

The  number  of  additions  and  multiplications  which  are  required  to 
implement  distance  and  other  methods  for  various  values  of  n  and  G,  are 
shown  in  Figures  4-1  and  4-2,  respectively.  From  Figures  4-1  and  4-2  it 
is  clear  that  if  n  is  greater  than  4,  the  distance  method  requires  less 
computation  when  compared  to  the  correlation  method.  In  order  to  have 
one  single  measure  of  comparison,  the  time  ratio  relating  multiplication 
and  addition  operations  can  be  used.  For  machines  like  the  IBM  370, 
multiplication  time  is  three  times  the  addition  time.  With  the  above 
assumption,  the  number  of  equivalent  additions  required  to  implement  cor¬ 
relation  and  distance  methods  for  various  values  of  n  and  G  is  shown  in 
Figure  4-3.  Figures  4-1,  4-2,  and  4-3  prove  that  the  distance  method  is 
computationally  more  efficient  than  the  correlation  method  if  the  number 
of  windows  is  large. 

Relation  between  Distance  Measures 
and  Moments 

Let  F  be  a  KxL  array  of  digital  pixels  which  may  assume  one  of  G 
possible  levels  on  the  gray  scale. 


Moments  method 


Number  of  Windows 

Figure  4-1.  Number  of  additions  to  register  n  windows. 


5.0x10 


Number  of  Windows 

Figure  4-2.  Number  of  multiplications  to  register  n  windows. 
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suoiqippy  }usi_eAinb3 


Figure  4-3.  Number  of  equivalent  additions  to  register  n  windows. 
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(4-56) 


1  <  F(X,Y)  <  G 

for  1<X<_K,l<YiL 

Let  S9  denote  the  set  of  pixels  with  pixel  value  equal  to 

the  number  of  elements  in  S^.  Now,  digital  image  F  can  be 

12  r 

the  union  of  sets  S  ,  S  ,  ....  S  . 

F  =  S1  U  S2  U  ...  U  SG 

Each  element,  which  is  just  a  point  in  the  two  dimensional 
the  set  S^,  is  uniquely  specified  by  its  coordinates. 


where 


and 


9  =  1.  2, 


,  for  i  =  1 ,  2, 


•  •  f 


G 


•  •  •  » 


n 


g 


The  set  can  be  considered  as  a  two  dimensional  discrete 

value  is  equal  to  g  at  P^,  p|,  ....  P^  and  zero  elsewhere 

q  a 

plane.  Therefore,  the  cpntral  moments  of  Sy  are  given  by 


g  and  n^  be 
considered  as 

(4-57) 

x-y  plane  of 


(4-58) 

function  whose 
in  the  x-y 

(4-59) 

(4-60) 


j 
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yg 


From  Equations  (4-59),  (4-60)  and  (4-61), 


(4-61 ) 


i°^02  .  I_  f  [(XJ_„,2  +  <*_*,*, 

“00  9  (=1 


(4-62) 


The  intraset  distance  of  the  set  S9  is  given  by  Equations  (4-36)  through 
(4-38). 


2  (T 

V>  w 


[(xf-x9)2  +  (vf-vS)2]) 


(4-63) 


Therefore, 


q  Q 

/ti20+y02y 

V  g  1 
u00 


(4-64) 


Equation  (4-64)  shows  the  relation  between  the  intraset  distance  and  mo¬ 
ments  of  the  set  S9.  Since  (u20+lJ02^  1S  1rwanant  under  translation  and 
rotation  [10],  the  intraset  distance  is  also  invariant  under  translation 
and  rotation.  If  each  point  P?  is  considered  as  a  point  mass  of  value 


9> 


q  a 
^20*^02  • 

,9 

u00 


the  square  of  the  radius  of  gyration  of  the  mass  distribu¬ 


tion  about  its  centroid.  Therefore, 


(4-65) 


where,  r  is  the  radius  of  gyration  of  S9  about  its  centroid. 
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The  interset  distance  between  the  sets  S9  and  Sh  is  defined  as  the 
distance  between  their  centroids. 

02(S9,Sh)  =  (X^-xV  +  (Y^-yV  (4-66) 

.  *10  Hll0^2  ,m01  m01  \2 

“  '  g  '  h  ;  '"g 

m00  m00  m00  m00 

Equation  (4-66)  shows  the  relation  between  the  interset  distance  between 
the  sets  S9  and  S*1,  and  their  moments.  Simulation  results  given  in  Refer¬ 
ence  [16]  show  that  matching  homogeneous  regions  separately  and  combining 
the  partial  results  additively  yields  sharper  correlation  peaks.  If  the 
image  is  composed  of  homogeneous  regions  as  described  in  Reference  [16], 
all  pixels  of  the  homogeneous  region  normally  fall  into  the  same  set  when 
the  image  is  segmented  based  on  pixel  values.  Computing  intraset  distance, 
in  some  sense,  is  the  same  as  processing  each  homogeneous  region  sepa¬ 
rately.  The  relative  location  of  homogeneous  regions  with  respect  to 
each  other  is  determined  by  interset  distances.  Therefore,  the  distance 
method  is  expected  to  perform  better  than  the  moments  method  if  the  scene 
is  composed  of  homogeneous  regions. 


Correlation  of  Adjacent  Pixels 
for  Image  Registration 


Let  W  be  a  digital  image  of  size  KxL  whose  pixels  can  assume  one 
of  G  possible  levels  on  the  gray  scale. 


for 


0  <  W(£,m)  <  G-l 
1  <  i  <  K  and  1  <  m  <  L 


(4-67) 


The  normalized  correlation  between  adjacent  pixels  of  the  izn  row 


of  the  digital  image  W  is  given  by  the  ensemble  average 
_  E{W(£,m)W(£,m+l )} 

p  -  2 

E{Vr(i,m)} 


(4-68) 


Since  there  are  L  elements  in  any  row  of  W,  p  can  be  approximated  by 
the  spatial  average 


t-U-  £  W(£,m)W(£,m+l) 
L  1  m-1 _ 

f  I  W2(£,m) 

L  m=l 


(4-69) 


where  P)l  always  lies  between  zero  and  one.  When  all  pixels  in  the  £th 
row  have  the  same  value,  p  is  one.  The  value  of  p  ,  in  some  sense,  is 

X-  JJ, 

related  to  the  difference  in  pixel  values  of  the  adjacent  pixels.  Values 
of  adjacent  pixels  are  highly  correlated  for  most  images  except  at  edge 
pixels.  If  pixel  W(£,m)  is  of  a  certain  gray  level,  then  the  adjacent 
pixel  W(£,m+1)  along  the  scan  line  £,  is  likely  to  have  a  similar  value. 
This  property  of  p£  has  been  used  in  image  coding  and  transmission  in  the 
past  [27,  pp.  278-2 81].  A  new  method  of  accomplishing  image  registration 
using  normalized  correlation  of  adjacent  pixels  of  rows  and  columns  of 
digital  images  is  presented  next. 


Feature  Vector  for  a 
Digital  Image 

Since  in  Equation  (4-69)  L/(L-1)  is  a  constant  and  can  be  dropped 
without  losing  any  information,  the  normalized  correlation  between  adja¬ 
cent  pixels  of  the  i**1  row  is  given  by 
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p*  '  L 


L-l 

X  W(*,m)W(i,,m+l) 
m=l 


I  VTU.m) 
m=l 


for 


l  =*1 ,  2,  ...»  K 


(4-70) 


Let  <rm  be  the  normalized  correlation  between  the  adjacent  pixels  of  the 
m^  column. 


m 


X  W(t,m)W(i+l ,m) 
z=l 
K  , 

X  knU.m) 

z=l 


(4-71) 


for  m  *  1,  2,  ....  L 

The  (K+L)-dimensional  feature  vector,  V  ,  of  the  digital  image  W  is  giv 

w 

en  by 


1  p2 


PK  °1  a2 


aL] 


(4-72) 


Multiple  Image  Registration 

A  method  of  accomplishing  multiple  image  registration  using  fea¬ 
ture  vectors  of  windows  and  subimages  of  the  search  area  S  is  described 
below. 

Step  1 :  Compute  the  feature  vector  for  each  of  the  n  windows  and  each 
of  (M-K+l )(N-L+1 )  subimages  of  S.  Let  and  denote  feature  vec¬ 
tors  of  the  window  W.  and  subimage  S_.  respectively. 

Step  2:  To  register  W,.  within  S,  find  the  subimaqe  S..*  .•*  whose  feature 

VJk 

vector  best  matches  that  of  W.  .  In  other  words,  find  S.*  •*  such  that 

VJk 
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IIV^v’Mcll2  <  I  |Vk-V1,j!  I2  (4-73) 

for  1  i  i  1  M-K+l,  i^i£,  and  1  <_  j  £  N-L+1 , 

k  i  i  2 

where  ||V  -V  j j  is  the  square  of  the  distance  between  the  vectors 
Vk  and  V1,j. 

To  register  windows  W^}  W^,  ....  Wn  within  S,  the  feature  vector 
for  n  windows  and  (M-K+l ) (N-L+1 )  subimages  are  computed  only  once  (Step 
1)  and  Step  2  is  repeated  n  times,  once  for  each  window.  Since  the  com¬ 
putation  associated  with  Step  2  is  small  compared  to  the  computation 
associated  with  Step  1,  this  method  is  more  efficient  for  multiple  image 
registration  than  the  correlation  and  SSDA  methods. 

Derivation  of  Multiplication 
and  Addition  Requirements 

Equations  for  the  number  of  multiplications  and  additions  in  terms 
of  M,  N,  K,  L  and  n  are  derived  in  the  seven  steps  below. 

1.  From  Equation  (4-70),  it  is  clear  that  2L  multipli¬ 
cations  and  2L-3  additions  must  be  performed  to  compute 
psl.  Similarly,  2K  multipl ications  and  2K-3  additions 
are  needed  to  compute  o  .  Since  there  are  K  rows  and 

L  columns,  a  total  of  4KL  multiplications  and  (4KL-3K- 
3L)  additions  are  required  to  compute  the  feature  vec¬ 
tor  for  a  KxL  array. 

2.  Let  p^i.j)  and  o  (i,j)  denote  the  normalized  corre¬ 
lation  between  adjacent  pixels  of  the  row  and  the 
m^  column  of  the  subimage  S.  .,  respectively.  In 

I  »J 
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order  to  compute  the  feature  vector  for  the  subimage 

S-j  1 1  4KL  multiplications  and  4KL-3K-3L  additions  are 

required  as  outlined  in  Step  1. 

11  12 

Once  V  *  is  computed,  V  *  can  be  computed  with 
very  few  arithmetic  operations  as  shown  below. 


L-l 


m 


L- 


m 


i  ) S-i  i(t,m+l) 

i=l  1,1  1 » 1 


l ^ 


U.m) 


At(l,l) 

'B^TTiy 


where 


A,  (1,1)  *  J  S,  nU.nOS,  ,(t,m+l) 
1  m=l  ’  1,1 

BjO.D  =  I  s*  ,(»,».) 

m=l  ’ 


Now, 


(4- 


(4- 

(4- 


1.-1 

PtO  *2)  =  ^ - - - - -  ( 

Z  S? 

m=l 

V1’1)‘Sl.l(t’1)Sl.l(t>2)  +  SK2(&,L-1)S1  2(&,l) 
~  Bjl.D-S^U.l)  +  sf^U.L) 

Frcwn  Equation  (4-77),  p^(l,2)  can  be  computed  from 
Pt(l,l)  by  performing  only  5  multiplications  and  4 
additions.  In  general,  p^i.j+l)  can  be  computed 
from  ot(1,j)  by  performing  5  multiplications  and  4 


additions. 
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Since  there  are  K  rows,  computation  of  p^i.j+l)  from 
p  (i,j)  for  i  -  1,  2,  ....  K  requires  a  total  of  5K 

Xi 

multiplications  and  4K  additions.  Also, 

am(i,j+l)  =  om+1 (i,j)  for  m  =  1,  2,  ....  L-l 

Therefore,  it  is  necessary  to  compute  o^(i,j+l)  which 
needs  2K  multiplications  and  2K-3  additions.  In  other 
words,  computation  of  the  feature  vector  from  V 

requires  7K  multiplications  and  6K-3  additions.  How¬ 
ever,  in  order  to  compute  V1*^  from  Vi,J,  Vi  and 

the  first  column  of  S.  .  must  be  stored.  Since  V1,J 

T  »  J 

is  a  (K+L)-dimensional  vector  and  S.  .  is  a  KxL  array, 

*  »J 

2K.+L  additional  memory  locations  are  required. 

4.  Similarly,  computation  of  Vi+^’J  from  V1 requires 
7L  multiplications  and  6L-3  additions  and  K+2L  addi¬ 
tional  memory  locations. 

5.  There  are  (M-K+1 ) (N-L+l )  allowable  reference  points 

in  the  search  area.  The  feature  vector  associated 
with  the  first  reference  point  (1,1)  is  computed  as 
described  in  Step  2.  This  requires  4KL  multipli¬ 
cations  and  4KL-3K-3L  additions.  ...  ,  VM~K+^ ’ 

are  computed  using  the  procedure  similar  to  the  one 
in  Step  3.  This  requires  (M-K)(6L-3)  additions  and 
7L(M-K)  multiplications.  For  line  i,  the  feature 
vectors  V^’^,  v^,  ...,  can  be  computed 

using  the  procedure  outlined  in  Step  3  with  7K(N-l) 
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multiplications  and  (6K-3)(N-L)  additions,  and  there 
are  (M-K+l)  such  lines.  Therefore,  to  compute  (M-K+l) 

(N-L+l )  feature  vectors: 

Number  of  multiplications  =  (4-79) 

4KL  +  7(M-K)[K(N-L)+L] 

Number  of  additions  =  (4-80) 

4KL  -  3K  -  3L  +  (M-K) (N-L) (6K-3)  +  (M-K) (6L-3) 

6.  The  computation  of  n  feature  vectors  for  n  windows 
needs  4KLn  multiplications  and  (4KL-3K-3L)n  additions. 

7.  (2K+2L-1 )  additions  and  (K+L)  multiplications  must  be 
performed  to  compute  ||Vk  -  V 1 ’ J | | 2 .  Since  there  are 
(M-K+l ) (N-L+l )  reference  points  and  n  windows, 

Number  of  multiplications  =  (K+L) (M-K+l ) (N-L+l )n  (4-81) 

Number  of  additions  =  (2K+2L-1 ) (M-K+l ) (N-L+l )n  (4-82) 

From  Step  5  through  Step  7, 

Total  number  of  multiplications  =  (4-83) 

(4KL+(M-K+1 ) (N-L+l ) (K+L) }n  +  4KL+7(M-K) [K(N-L)+L] 

Total  number  of  additions  =  (4-84) 

{(2K+2L-1 ) (M-K+l ) (N-L+l )+(4KL-3K-3L) }n  + 

(M-K) (N-L) (6K-3 )+(M-K) (6L-3 )+4KL-3K-3L 

For  M  =  240,  N  =  256,  K  =  32,  and  L  =  32, 

Total  number  of  multiplications  =  (4-85) 

301 3696n  +  10487296 


4 


Figure  4-4,  Number  of  additions  to  register  n  windows 


Figure  4-5.  Number  of  multiplications  to  register  n  windows. 


0l*0‘8 


Figure  4-6.  Number  of  equivalent  additions  to  register  n  windows. 
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Total  number  of  additions  ■  (4-86) 

597607 9n  +  8849104 

The  nionber  of  additions  and  multiplications  which  are  required  to 
Implement  several  registration  methods  for  various  values  of  n  are  shown 
In  Figures  (4-4)  and  (4-5),  it  Is  obvious  that  the  new  method  based  on 
correlation  between  adjacent  pixels  is  computationally  more  efficient 
than  any  of  the  methods  previously  considered.  Since  or  V*+1*^ 

can  be  computed  from  with  few  arithmetic  operations,  this  method 
Is  promising  for  real  time  Implementation.  A  few  simulations  were  run 
using  Images  from  similar  sensors  and  the  above  method  was  successful. 


Feature  Extraction  Technique  for 


strati  on 


When  two  Images  do  not  differ  In  pixel  resolution  and  rotation, 
the  method  most  widely  used  for  Image  registration  Is  cross-correlation. 
The  elements  of  the  normalized  cross-correlation  surface  are  defined  to 


K  l 


R(1»j) 


T  £  W(i,m)  S.  ,(i,m) 
m»l 


K  L 


Cl  I  W2(z,m)rtf  T  S2  (t.m)]1 
m-1  &  m*l  1,J 


(4-87) 


1  <  i  <  M-K+l ,  1  <  j  <  N-l+1 


In  general,  the  amount  of  computation  associated,  with  any  simi¬ 
larity  detection  method  Is  proportional  to  the  number  of  pixels  In  the 
window  and  the  number  of  pixels  In  the  al lowable  search  area.  For  the 
cross-correlation  method,  each  of  the  KL  pixels  In  W  Is  compared  with 


the  corresponding  KL  pixels  in  j  to  compute  R(i,j).  Since  the 
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correlation  function  has  (M-K+1 )(N-L+1 )  elements,  a  total  of  KL(M-K+1) 
(N-L+l)  pairs  of  pixels  are  compared.  Thus  the  computation  time  Is 
roughly  proportional  to  Kt(M-K+l) (N-L+l).  Schemes  such  as  "two-stage 
template  matching"  and  "course-fine  template  matching"  have  been  pro¬ 
posed  to  speedup  correlation  method  and  are  described  In  Chapter  I.  A 
feature  extraction  technique  presented  In  this  Section  selects  a  set  of 
pixels,  H1,  from  the  window  W  to  be  used  In  correlation  and  thus  achieves 
significant  savings  In  computation  with  little  effect  on  the  correlator 
accuracy  and  reliability. 

All  methods  which  speedup  correlation  accomplish  the  task  by 
somehow  reducing  the  total  number  of  pixel  pair  comparisons  during  com¬ 
putation  of  the  correlation  surface.  Feature  extraction,  one  of  the 
fundamental  methods  for  data  compression  In  the  field  of  pattern  recog¬ 
nition  can  be  used  to  accomplish  the  same.  An  Ideal  feature  Is  re¬ 
quired  to  have  the  following  properties: 

1.  The  feature  should  retain,  from  the  original  pattern, 
as  much  information  as  possible. 

2.  The  feature  should  accomplish  as  much  data  reduction 
as  possible. 

3.  The  feature  should  be  Invariant  or  depend  on  some  In¬ 
variant  properties  of  the  original  pattern  In  a  known 
way. 

In  practice  the  first  two  above  are  conflicting  properties.  However, 
striking  a  balance  with  consistent  and  acceptable  accuracy  for  recog¬ 
nition  of  the  original  pattern  can  be  accomplished. 


100 

Consider  the  subset  of  pixels,  W,  from  W  as  a  feature  based  on 
the  mean  and  standard  deviation  of  pixel  values  of  W  as  shown  in  Figure 
4-7.  W  is  a  set  of  all  pixels  in  W  such  that  W(t,m)  is  either  greater 
than  the  mean  plus  a  constant  p  times  the  standard  deviation  (y+pcr)  or 
less  than  (y-po),  where  p  is  the  scale  factor  greater  than  zero.  This 
feature  retains  pixels  from  W  whose  pixel  values  are  relatively  low  or 
high  as  compared  to  the  mean.  Data  reduction  is  accomplished  by  delete- 
Ing  all  pixels  in  W  whose  pixel  values  range  from  (y-po)  to  (y+po).  The 
amount  of  information  retained  and  the  number  of  pixels  n  In  W,  depends 
on  the  scale  factor  p  and  therefore  is  controllable.  Finally,  the  fea¬ 
ture  depends  on  statistics  of  the  window  in  a  known  way. 

Elements  of  the  correlation  surface  are  computed  by  comparing 
only  the  n  pixels  from  W  which  belong  to  feature  set  W‘  with  their  cor¬ 
responding  pixels  In  the  subimages  of  S. 

Z  WU,m)  S^jU.m) 

R(1,J)  »  - - - g -  (4-88) 

[  I  l  S,  Ji, m)]*5 

(&)  U,m) 

for  all  (t,m)  such  that  W(t,m)  belongs  to  W  and  for 
1  <  i  <  M-K+l ,  1  <  j  <  N-L+l 

Now  the  computation  time  is  roughly  proportional  to  n(M-K+l) 
(N-L+l).  The  percent  savings  in  computation  time  depends  on  the  histo¬ 
gram  of  W  and  the  scale  factor  p  and  therefore  is  scene  dependent. 
However,  simulation  results  presented  In  the  next  section  indicate  a 
savings  of  50  to  75  percent  if  p  Is  set  to  one.  The  percent  savings  in 
computation  time,  T,  can  be  computed  using  the  relation 


T  *  (i^)*100* 
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(4-89) 


In  real  world  problems,  n  depends  on  the  scene  and  nothing  much 
can  be  said  about  It  without  knowledge  of  the  histogram  of  W.  However, 
In  order  to  gain  some  Insight  Into  this  method,  consider  the  following 
two  examples. 

Example  1:  Let  W  be  an  Image  whose  pixel  values  assume  a  unlmodal 
gausslan  distribution  with  mean  p  and  standard  deviation  o;  i.e.. 


2 

p(x) - -  exp[-kl£L]  (4-90) 

^2 *■  a  Za 

Then, 

ktPa 

number  of  pixels  in  H'  .  ,  / 

number  of  pixels  in  H  1  ~  /  !><x>dx  (4-9,) 

x*y^>a 

.  t  .  [erf(Hiep.)  .  erf(Hl£pL)] 

■1-2  erf(p) 


where 


erf(p) 


dx 


Therefore,  T  *  200  erf(p)% 


A  plot  of  T  versus  p  Is  shown  in  Figure  4-8.  For  small  values  of 
p,  the  curve  is  almost  linear  Indicating  rapid  reduction  in  computation 
time  as  p  is  Increased  from  zero.  For  large  values  of  p  the  curve 
becomes  flat  yielding  small  savings  In  computation  for  corresoonding 
Increases  In  p. 


0  0.5  1.0  1.5  2.0  2.5  3.0 

Scale  factor 

Figure  4-8.  Plot  of  percent  savings  In 

computation  T  versus  scale  factor 


103 


Example  2:  Let  W  be  a  blmodal  Image  which  means  that  the  histogram  of 
U  Is  a  mixture  of  two  unlmodal  histograms  p-j(x)  and  P2(x).  u-j  Is  the 
mean  of  p-j(x)  and  a.j  Is  the  standard  deviation  about  uj.  Similarly  u2 
Is  the  mean  of  p2(x)  and  a2  Is  the  standard  deviation  about  u2.  If 
and  P2  are  aprlorl  probabilities  of  two  principal  brightness  levels, 
then  the  histogram  of  W,  p(x),  can  be  expressed  as  [26,  pp.  326*328] 


p(x)  *  P1  p-,(x)  +  P2  p2(x) 
where  P^  +  P2  *  1 

For  the  gausslan  case 


Pt  (x) 


o-j 


exp[ - 


2 


P2(x) 


r‘(x'w2)  i 
exp[ - r—] 

2  at 


Let  v  and  a  be  the  mean  and  variance  of  pixel  values  of  W. 


(4-92) 


(4-93) 


(4-94) 


u  *  E[x]  *  Pjii'j  +  P2n2 


(4-95) 


E[x2] 


CP1P1 (x)  +  P2P2(x)]x2  dx 


(4-96) 


Pj  (w2  +  o2)  +  P2(u2  +  ff2) 


2 

The  variance  a  Is  given  by: 
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?  -  E[x2]  -  y2 

*  Pl^l  +  °l^  +  P2^w2  +  a2^  '  ^Plyl  +  P2U2^ 

*  Plul  "  Plwl  +  P2m2  “  P2W2  “  2P1P2w1u2  +  Plal  +  P2 
«  P}u20  -  P,)  +  P,p?(1  -  P9)  -  2?,?^,  +  PlCT?  + 


1'  -  r2M2 


lr2“r2  rl°l 


PlP2(l,l  *  w2)2  +  Plcl  +  P2a2 


The  percent  savings  in  computation  T  is  given  by: 


p+pa 

T  *  J  p(x)  dx 

p-pa 

p+pa 


/  CPtP^x)  +  P2P2(x)]  dx 


p-pa 


p+po 


u+pa 


*  P1  J  P-|(x)  dx  +  P2  J  p2(x) 


dx 


p-pc 


p-ap 


p+po-p.  y-pp-p, 

'PjlerfC — - — -]  -  erf [ — — — 4> 


p+OO-p, 

+  P9{erf[—  r---]  -  erft 
c  a0 


C1 

p-pa-p2 


]} 


Substituting  Equation  (4-95)  Into  the  above  equation  yields 

V 


Mpo-WiJ+PO  P^p^-pj-po 

T  -  P,{erf[  -  _  1 - ]  -  erf[— — \ — - - ■]} 


'1 


’1 


P,(p,-p?)+pc  P-, (p,-y9)-po 

+  P?{erf[— ! — ^ — = - ]  -  erf[-— '  2 - •]} 
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where  a2  ■  piP2^wi"w2^2  +  Plffl  +  p2a2 

Notice  T  depends  on  the  parameters  y-j,  a]  and  of  the  image  and  of 
course  on  the  scale  factor  p.  T  is  scene  dependent  and  general  conclu¬ 
sions  cannot  be  drawn  without  knowledge  of  the  histogram  of  W. 


Simulation  Results 

Scenes  used  for  simulation  were  obtained  from  sensors  sensitive  in 
the  visual  spectrum  (day  TV  sensors).  Because  of  the  difference  in  the 
sensors,  the  two  Images  were  preprocessed  such  that  they  have  the  same 
spatial  resolution.  An  algorithm  to  accomplish  this  is  given  in  Refer¬ 
ence  [4].  W  is  a  32x32  array  of  pixels  and  S  Is  a  120x120  array  of 
pixels  extracted  from  preprocessed  high  and  low  resolution  images,  re¬ 
spectively.  Since  a  one-bit  or  two-level  correlator  is  used  in  this 
work,  U'  and  S  must  be  quantized  to  two  levels.  The  mean,  u,  and  stan¬ 
dard  deviation,  o,  of  W  are  first  computed.  Pixels  are  quantized  to  1 
or  0  if  they  are  above  y+pc  or  below  y-pa,  respectively.  The  locations 
of  pixels  within  U  with  values  between  y-po  and  y+pa  are  stored  and  mask¬ 
ed  out  of  the  correlation  process.  Also,  each  subimage,  Si  ^  of  S  was 
quanitzed  to  ones  or  zeros  depending  on  whether  the  pixel  value  was 
above  or  below  the  average  value  of  the  subimage,  y.  ,,  as  given  below 

•  *  J 


1 

Sj  .,U,m)  *  { 

•  »j  0 


,  If  S^jU.m)  >  y1fj 
,  otherwise 


(4-99) 


A  number  of  simulations  were  run  using  typical  scenes  for  different 
values  of  p.  As  a  figure  of  merit,  the  correlation  surface  slgnal-to- 
noise  ratio  was  computed  as 
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SNR  *  Ry  ~-av3  ,  (4-100) 

sigma 

where  is  the  maximum  value  of  correlation  surface, 

Ravg  avera9e  value  correlation  surface, 

^sigma  stan(^ar(^  deviation  of  correlation  surface. 

Correlation  was  found  successful  with  the  highest  value  in  the 
correlation  surface  Indicating  the  true  registration  for  values  of  p 
ranging  from  0  to  1.4.  When  p  was  Increased  beyond  1.4,  correlation 
was  unsuccessful.  Percent  savings  in  computation  time  and  signal -to- 
noise  ratio  of  the  correlation  surface  for  various  values  of  p  ranging 
from  zero  to  one  are  tabulated  in  Table  4-1  and  Table  4-2,  respectively. 
From  the  above  simulation  results,  it  is  concluded  that  substantial 
savings  in  computation  (50  to  75  percent)  is  accomplished  with  little 
effect  on  correlation  accuracy  and  reliability. 

Generalization 

Since  the  new  feature  extraction  technique  reduces  the  amount  of 
data  to  be  processed  without  altering  the  structure,  it  is  applicable 
not  only  to  correlation  but  to  many  other  similarity  detection  methods. 
The  following  generalizations  can  be  made: 

1.  By  using  n  elements  of  W'  as  test  points  in  sequential 
similarity  detection  algorithm  suggested  by  Barnea 
and  Silverman,  the  problem  of  digital  image  regis¬ 
tration  reduces  to  the  problem  of  finding  Siit  .*  such 

■  » J 

that 
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Table  4-1 

Percent  savings  in  computation. 


Scene 

p-0.00 

P-0.25 

p-0. 50 

p-0.75 

p-1.00 

1 

00.00 

16.89 

33.30 

57.92 

75.78 

2 

00.00 

13.21 

28.90 

48.82 

71.38 

3 

00.00 

3.61 

12.03 

27.63 

51.30 

4 

00.00 

19.02 

27.52 

53.21 

70.50 

5 

00.00 

20.00 

37.42 

55.34 

73.24 

Table  4-2 

SNR  of 

correlation  surface. 

Scene 

p=0. 00 

p-0.25 

p- 0.50 

p=0.75 

p-1.00 

1 

5.885 

5.609 

4.907 

4.669 

3.934 

2 

7.175 

6.101 

5.045 

4.356 

4.321 

3 

6.775 

6.769 

6.650 

6.347 

5.236 

4 

7.619 

7.20 

6.64 

6.13 

5.31 

5 

8.197 

8.14 

7.91 

7.51 

6.85 
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(4-101) 


In  (4-101) 


2. 


3. 


e.*  ■?*  <  e<  -s  f°r  1  1  1  <  M-K+l  and  i  t  i* 

•  Id  * 

1  <  j  <  N-L+l  and  j  t  j* 


e.  •  can  be  computed  using  Equation  (4-102) 

■  *J 


I  |S.  .(i,n)  -  S.  .  -  W(t,m)  +  W| 
(Mi)  1,J  1,J 


where 


=  1  I 

n  / 


(t,m) 


S.  .(*, m) 
* 


and 

W  -  i  I  M(£,m) 

"  (A) 

for  all  (t,m)  such  that  W(t,m)  belongs  to  W'  and  for 
1  <  i  <  M-K+l ,  1  <  j  <  N-L+l 

A  few  simulations  were  run  and  the  above  method  was 
found  successful . 

For  the  moments  method,  moments  can  be  computed  based 
on  n  elements  of  W*  rather  than  all  the  KL  elements 
of  W. 

This  method  can  be  combined  with  the  improved  method 
for  correlating  similar  sensor  images,  suggested  in 
Reference  [30]  to  improve  the  probability  of  finding 
the  true  peak. 

Step  1 ;  Compute  the  Cross-Correlation  surface  by 
comparing  only  the  n  pixels  from  W  which  belong  to 
the  feature  set  W  with  their  corresponding  pixels 
in  the  subimages  of  S  (using  Equation  (4-38)). 


(4-102) 


(4-103) 


(4-104) 
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Step  2:  Identify  a  predetermined  number  of  highest 
peaks  and  coordinates  of  their  occurence  from  the 

correlation  surface.  Let  (^2**^2^  *’*  * 

(I,  ..j.  .)  and  be  the  coordinates  of  the  first 

'  lc-i  k-l  K  K 

k  peaks.  It  is  assumed  that  the  true  registration 
point  is  one  among  them.  Recompute  R(  1^ , Ji ) •  Rflg*^* 

...  ,  and  a11  the 

KL  pixels  of  W  with  the  corresponding  pixels  in  the 
subimages  of  S  beginning  at  (I-|»Jj)*  »  •••  * 

(Ik-rJk-l)  and  ^k’V*  resPect1‘vely  (usin9  Equation 
(4-87)).  Savings  in  computation  is  accomplished  in 
Step  1.  Step  2  increases  the  probability  of  finding 
the  true  peak  and  reduces  the  possibility  of  false 
registration. 


HOCKXUB  PASS  BUHL-NO*  FILMED 


V.  CONCLUSIONS  AND  RECOMMENDATIONS 

Most  of  the  conclusions  and  recomnendations  presented  In  this 
chapter  have  been  given  and  justified  within  the  first  five  chapters 
of  this  dissertation. 


Conclusions 

1.  Various  methods  for  accomplishing  digital  image  registration  are 
presented  in  Chapter  II.  The  most  commonly  used  method  is  cross¬ 
correlation.  There  exist  two  independent  ways  of  computing  the 
correlation  surface  (Direct  method  and  the  FFT  method).  Even 
though  the  FFT  method  requires  less  computation  compared  to  the 
direct  method,  it  requires  a  large  amount  of  memory  for  software 
Implementation  and  Is  very  complex  for  real  time  hardware  imple¬ 
mentation.  The  direct  method  requires  less  memory  compared  to  the 
FFT  method,  involves  no  complex  multiplications  or  complex  addi¬ 
tions  and  can  be  easily  Implemented  using  digital  hardware. 

2.  Even  though  the  vector  correlation  is  expected  to  yield  better 
performance  than  the  standard  correlation  algorithm,  its  use  is 
limited  by  the  large  amount  of  computation  required  to  implement 
the  method  (more  than  twice  the  computation  required  by  the  stan¬ 
dard  correlation  algorithm). 

3.  Due  to  the  difficulty  encountered  in  the  automatic  segmentation  of 
digital  images  into  homogeneous  regions,  feature  matching 
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correlation  and  hybrid  correlation  algorithms  may  not  be  of  any 
significant  use  (especially  for  hardware  implementation). 

4.  Sequential  similarity  detection  algorithm  computes  error  surface 
as  a  measure  of  dissimilarity  between  the  window  and  subimages  of 
the  search  area.  Since  this  method  requires  only  addition  and  a 
few  division  operations.  Implementation  is  simpler  than  that  of 
the  correlation  method  (addition  or  subtraction. is  simpler  than 
multiplication  or  division). 

5.  All  the  three  popular  techniques  which  speedup  template  matching 
(two-stage  temolate  matching,  course-fine  temolate  matching  and 
hierarchical  search  method),  accomplish  savings  in  computation 

by  reducing  the  total  number  of  pixel  pair  comparisons.  The  meth¬ 
ods  are  computationally  more  efficient  in  terms  of  the  number  of 
arithmetic  operations  required  (software  Implementation)  but  may 
not  enjoy  any  advantage  in  real-time  implementation  using  special 
purpose  hardware. 

6.  The  necessity  of  a  highly  efficient  algorithm  to  transform  digital 
Images  to  binary  Images  (edge  and  non-edge  pixels)  and  a  connec¬ 
tivity  test  to  Identify  true  straight  line  edges  (composed  of  adja¬ 
cent  pixels)  makes  the  use  of  the  Hough  transformation  for  digital 
Image  registration  less  attractive. 

7.  The  moments  of  an  Image  or  subimage  can  be  easily  computed  by  per¬ 
forming  multiplications  and  additions  only.  When  the  window  and 
the  search  area  do  not  differ  In  rotation  or  pixel  spatial  resolu¬ 
tion,  It  Is  not  necessary  to  compute  Hu's  invariant  moments  and  the 
ordinary  moment  sequence  can  be  directly  used  for  scene  matching. 
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8.  For  the  reasons  given  In  1  through  7,  It  Is  concluded  that  the 
standard  correlation  algorithm,  the  SSOA  and  the  moments  method 
are  more  promising  using  present  state-of-the-art  hardware. 

9.  In  Chapter  III,  the  computational  efficiency  of  the  correlation, 
the  SSDA  and  the  moments  methods  Is  compared  from  software  as  well 
as  hardware  points  of  view,  Independently,  for  the  multiple  Image 
registration  problem.  It  Is  found  that  the  moments  method  takes 
less  computation  time  If  Implemented  In  software  and  less  hardware 
for  real-time  Implementation  If  the  number  of  windows  Is  sufficient¬ 
ly  large.  Moments  of  any  order  can  be  computed  with  one  level  of 
multiplications  and  one  level  of  additions.  Therefore,  the  moments 
method  Is  promising  for  real-time  Implementation. 

10.  In  order  to  obtain  better  results  from  any  of  the  Image  registration 
algorithms,  the  window  and  each  subimage  of  the  search  area  should 
be  preprocessed  to  have  zero  mean  and  unity  standard  deviation. 

This  Is  called  Intensity  level  normalization  and  would  require  too 
much  additional  computation.  In  Chapter  IV,  It  is  shown  that  In¬ 
tensity  level  normalization  can  be  embedded  within  the  moments 
method  with  almost  no  additional  computation.  This  makes  the  mo¬ 
ments  method  more  accurate  without  excessive  additional  computation. 

11.  Two  new  feature  matching  algorithms,  one  based  on  Interset  and  In¬ 
traset  distances,  and  the  other  based  on  correlation  of  adjacent 
pixels  are  developed  In  Chapter  IV.  The  distance  method  Is  com¬ 
putationally  more  efficient  (In  terms  of  the  number  of  arithmetic 
operations)  than  the  correlation  or  the  SSDA  method  If  the  number 
of  windows  Is  greater  than  three. 
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12.  It  his  bun  shown  that  the  Intrasct  distinct  of  a  set  Is  a  function 
of  zero  and  second  order  central  moments  of  the  set.  The  Inter set 
distance  between  two  sets  Is  a  function  of  zero  and  first  order  mo¬ 
ments  of  both  the  sets.  When  an  Image  Is  composed  of  homogeneous 
regions,  all  pixels  of  a  homogeneous  region  normally  fall  Into  the 
same  set  when  the  Image  Is  segmented  based  on  pixel  values.  Com* 
putlng  Intraset  and  Interset  distances.  In  some  sense.  Is  the  same 
as  computing  moments  of  each  homogeneous  region,  separately.  In 
general,  processing  each  homogeneous  region  separately  and  combin¬ 
ing  the  partial  results  addltlvely  yields  better  performance. 
Therefore,  the  distance  method  Is  expected  to  perform  better  than 
the  moments  method  If  the  scene  Is  composed  of  homogeneous  regions. 

13.  For  multiple  Image  registration,  the  algorithm  based  on  correlation 
of  adjacent  pixels  Is  computationally  more  efficient  than  any  of 
the  algorithms  considered.  Even  for  single  Image  registration, 
this  method  requires  less  number  of  arithmetic  operations  than  the 
standard  correlation  algorithm.  Since  the  feature  vector  associated 
with  reference  point  (1,j+l)  or  (1+1 ,j)  can  be  computed  from  the 
feature  vector  associated  with  the  reference  point  (1,j)  with  very 
few  arithmetic  operations,  this  method  Is  very  promising  for  real¬ 
time  Implementation. 

14.  From  9,  11  and  13,  It  Is  concluded  that  an  algorithm  which  Is  com¬ 
putationally  efficient  for  single  Image  registration  may  not  be 
efficient  for  multiple  Image  registration.  It  Is  also  concluded 
that  the  computation  for  template  matching  algorithms  Is  directly 
proportional  to  the  number  of  windows  whereas  In  feature  matching 
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algorithms  features  are  extracted  for  all  sublnages  and  windows 
only  once,  and  the  Hatching  procedure  Is  repeated  once  for  each 
window.  Since  the  computation  required  to  natch  the  features  Is 
negligible  compared  to  that  required  to  compute  features.  In  general, 
feature  matching  algorlthas  are  expected  to  be  more  efficient  than 
tmqriate  matching  algorithms  If  the  number  of  windows  Is  sufficien¬ 
tly  Urge. 

15.  There  are  many  other  feature  matching  Image  registration  algorithms 
reported  In  the  literature  that  have  not  been  mentioned  in  this  re¬ 
port.  The  reason  for  eliminating  these  algorithms  was  either  they 
were  too  computationally  complex  to  be  Implemented  In  real-time 
hardware  In  the  near  future  or  they  were  not  as  computationally  ac¬ 
curate  as  the  methods  presented.  Special  purpose  VLSIC  developed 
at  some  future  date  may  make  some  of  these  algorithms  feasible. 

16.  A  feature  extraction  technique  based  on  the  mean  and  standard  de¬ 
viation  of  pixel  values  of  the  window  accomplishes  50  to  75  percent 
savings  In  computation  with  very  little  effect  on  registration  ac¬ 
curacy.  Since  the  feature  extraction  technique  reduces  the  amount 
of  data  to  be  processed  without  altering  the  structure.  It  Is  ap¬ 
plicable  not  only  to  correlation  but  to  many  other  Image  regis¬ 
tration  algorithms. 

Recoamendatlons  for  Future  Work 

1.  Even  though  the  correlation,  the  SSDA  and  the  moments  methods  are 
more  promising  using  present  state-of-the-art  hardware,  very  large 
scale  Integrated  circuits  developed  to  perform  a  specialized  task 
may  at  some  future  date  make  any  of  the  methods  discussed  In  Chapter 
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II  fusible.  Because  of  this  possibility,  a  follow-on  program  should 
Investigate  the  computational  accuracy  of  these  methods. 

2.  It  may  be  possible  to  obtain  good  registration  results  using  very 

few  moments.  This  trade-off  should  be  Investigated  through  Simula- 

0 

tlon  of  typical  military-type  digitized  scenes. 

3.  If  the  Images  are  quantized  to  two  levels  (0  and  1)  or  three  levels 
(-1,  0  and  +1),  the  computation  of  moments  Involves  no  multiplica¬ 
tions.  Therefore,  performance  of  the  moments  method  with  two  and 
three  level  Images  as  Inputs  should  be  studied. 

4.  Accuracy,  reliability  and  sensitivity  to  noise  of  all  three  feature 
matching  methods  In  Chapter  IV  should  be  determined  through  simula¬ 
tion. 

5.  The  effect  of  quantization  of  Images  to  two  or  three  levels  on  the 
Performance  of  feature  matching  algorithms  should  be  Investigated. 

6.  The  potential  of  the  method  based  on  correlation  of  adjacent  pixels 
for  real-time  hardware  Implementation  should  be  studied. 
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